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The  nonlinear  acoustic  propagation  effects  require  a 
numerical  solution  in  the  time  domain.  To  model  a 
porous  ground  surface,  which  in  the  frequency  domain 
exhibits  a  finite  impedance,  the  linear  phenomenological 
porous  model  of  Morse  and  Ingard  was  used.  The 
phenomenological  equations  were  solved  in  the  time 
domain  for  coupling  with  the  time  domain  propagation 
solution  in  the  air. 

The  numerical  solution  is  found  through  the  method  of 
finite  differences.  Two  kinds  of  numerical  absorbing 
boundary  conditions  were  developed  for  the  air  propaga¬ 
tion  equations  to  truncate  the  physical  domain  for  solution 
on  a  computer.  Radiation  conditions  first  were  used  on 
those  sides  of  the  domain  where  there  were  outgoing 
waves.  Characteristic  boundary  conditions  are  employed 
near  the  acoustic  surface.  Curves  of  pressure  amplifi¬ 
cation  versus  incident  angle  for  waves  obliquely  incident 
on  the  hard  and  porous  surfaces  were  produced. 

The  model  predicted  that  near  grazing  finite  amplitude 
acoustic  blast  waves  decay  with  distance  over  hard 
surfaces  as  r ' 2.  For  propagation  over  porous  ground  sur¬ 
face,  the  model  predicted  that  this  surface  decreased  the 
decay  rate  with  distance  for  the  larger  blasts  compared  to 
the  rate  expected  in  the  linear  acoustics  limit. 
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1.  INTRODUCTION 

What  do  the  following  situations  have  in  common:  a  wailing  civil  defense  siren, 
a  refrigerating  thermoacoustic  engine,  a  jet  airplane  taking  off  from  an  airport,  a 
hyperthermia  treatment  for  a  cancerous  tumor,  an  outdoor  blast  from  a  small  explo¬ 
sion,  a  lithotripsy  treatment  for  breaking  up  a  kidney  stone,  a  parametric  array  on  a 
submarine  sending  a  low  frequency  sonar  signal,  and  a  4th  of  July  fireworks  display? 
In  each  of  these  circumstances  sound  waves  of  high  amplitude  are  present.  The 
civil  defense  siren,  the  jet  taking  off,  the  small  explosion,  and  the  4th  of  July 
fireworks  are  all  loud,  very  loud  if  you  are  close  to  them.  You  may  not  be  able  to 
hear  directly  with  your  ears  a  thermo  acoustic  engine,  or  a  hyperthermia  treatment, 
or  the  sound  of  a  parametric  array.  But  each  of  these  cases  is  analogously  loud, 
involving  the  same  physical  processes  as  the  loud  sound  waves  outdoors. 

This  dissertation  presents  a  method  for  calculating  numerically  how  the  sound 
from  loud  sources  propagates  through  an  acoustic  medium.  The  specific  problem 
addressed  in  this  work  is  the  propagation  of  blast  noise  in  a  homogeneous  atmo¬ 
sphere  near  hard  and  porous  ground  surfaces.  To  verify  the  computer  model 
developed  for  blast  sounds,  the  sound  of  an  electric  spark  pulse  propagating  in  the 
free  field  and  near  hard  surfaces  is  also  examined  in  detail.  The  spark  problem 
occurs  on  a  scale  measured  in  centimeters,  and  the  blast  noise  occurs  over  distances 
measured  in  many  meters,  but  again  the  physical  problems  in  the  two  situations  are 
nearly  identical.  The  only  difference  between  the  two  is  that  the  amount  of  attenua¬ 
tion  in  the  air  is  greater  for  sparks. 
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The  measurement  and  characterization  of  the  impulse  noise  from  blasts  are 
critical  to  the  prediction  of  the  environmental  impacts  of  military  operations.  One 
must  take  measurements  close  enough  to  a  source  to  eliminate  metecro logical  varia¬ 
tions  but  far  enough  away  that  finite  amplitude  wave  effects  do  not  dominate.  In 
most  cases  this  restriction  implies  that  one  must  make  these  measurements  over 
natural  outdoor  surfaces.  The  effects  of  the  finite  impedance  of  the  ground  on  the 
propagation  of  linear  continuous  waves  are  profound,  and  this  fact  is  widely 
accepted. 

The  techniques  for  predicting  sound  propagation  out'oors,  using  the 
infinitesimal  pressure  amplitude  assumption,  linear  acoustic  methods,  are  well  esta¬ 
blished.  The  fast  field  program  (the  FFP)  and  the  parabolic  equation  method  (the 
PE)  are  currently  the  most  prevalent  computational  approaches  using  the  linear 
theory.1-6  These  methods,  however,  do  not  accurately  model  the  physics  when  the 
sounds  become  very  loud,  over  150  dB  referenced  to  20  //Pa,  and  the  amplitudes  of 
the  pressure  variations  making  up  the  sound  become  finite  instead  of  infinitesimal. 

This  dissertation  develops  a  numerical  model  to  predict  sound  propagation  at 
high  acoustic  amplitudes  where  the  infinitesimal  theory  breaks  down.  The  model 
will  include  the  contributions  of  geometrical  spreading,  ground  effects,  attenuation, 
and  nonlinear  distortion  to  blast  pulse  propagation.  With  this  model  it  will  be  possi¬ 
ble  to  make  predictions  for  how  far  a  receiver  should  be  from  a  blast  of  a  certain 
charge  size  for  the  infinitesimal  acoustic  propagation  theory  to  be  considered  valid. 
This  knowledge  will  allow  measurements  to  be  made  as  close  to  blast  sources  as 
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possible  to  minimize  meteorological  effects  but  far  enough  away  that  nonlinear 
acoustical  effects  will  not  dominate. 

The  remainder  of  this  thesis  is  divided  into  chapters.  The  next  five  chapters 
will  present  the  background  for  the  physical  problem  and  detail  the  equations  which 
describe  it.  Chapter  2  will  review  some  of  the  other  methods  researchers  are  using 
to  explore  nonlinear  acoustics  and  some  of  the  fundamentals  governing  finite  ampli¬ 
tude  waves  assuming  some  knowledge  of  linear  acoustic  theory.  Chapter  3  will 
review  some  of  the  considerations  involved  in  outdoor  sound  propagation,  especially 
the  effect  of  ground  impedance.  Chapter  4  will  detail  the  finite  amplitude  equations 
for  a  fluid  such  as  air  and  will  manipulate  these  equations  into  a  form  suitable  for  a 
numerical  solution.  Chapter  5  specifies  a  low  frequency  phenomenological  model 
for  a  porous  ground  surface.  The  initial  conditions  appropriate  for  blasts  and  spark 
pulses  are  described  in  Chapter  6. 

Chapters  7  through  9  give  the  numerical  methods  used  to  simulate  the  physical 
problem.  Chapter  7  presents  the  finite  difference  procedures  employed  for  the  air 
propagation  numerical  solution  and  some  motivation  as  to  why  the  method  was 
chosen.  Chapter  8  describes  the  numerical  absorbing  boundary  conditions  used  to 
truncate  the  physical  domain  of  the  atmosphere  for  simulation  on  a  computer  with 
finite  memory.  The  numerical  solution  found  for  the  porous  medium  ground  sur¬ 
face  is  in  Chapter  9. 

Chapter  10  gives  examples  to  show  that  the  numerical  simulation  method  is 
operating  correctly.  Here,  electric  spark  free  field  propagation,  normal  reflection, 
and  oblique  reflection  from  hard  surfaces  are  investigated  in  detail.  This  chapter 
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also  discusses  the  values  of  artificial  viscosity  needed  for  the  blast  wave  propagation 
runs. 

Chapter  11  presents  the  results  of  the  numerical  simulations  for  blast  waves 
propagating  over  the  hard  and  porous  medium  surfaces,  both  for  oblique  and  near 
grazing  geometries.  Chapter  12  concludes  the  main  body  of  this  thesis  and  points  to 
possible  future  work  using  and  extending  the  methods  given  here.  Appendix  A 
applies  finite  amplitude  acoustic  propagation  outdoors  as  addressed  in  this  disserta¬ 
tion  to  music,  specifically  to  performances  of  Tchaikovsky’s  Overture  1812  with  can¬ 
nons.  Appendix  B  describes  the  computer  program  mcnalp  developed  for  the 
numerical  simulations. 
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2.  AN  INTRODUCTION  TO  NONLINEAR  ACOUSTICS 

The  theory'  of  finite  amplitude  sound  waves,  or  nonlinear  acoustics,  is  a  well 
understood  subject.  Several  recent  books  have  been  written  on  the  topic,  or  have 
chapters  devoted  to  it  as  well  as  summary  articles  in  journals.'-10  These  references 
can  provide  a  history  of  the  subject,  which  will  no .  be  detailed  here.  The  following 
paragraphs  review  some  of  the  results  of  this  theory. 

Sound  waves  are  dilatational,  or  longitudinal,  waves  which  propagate  in 
compressible  fluids  and  solids.  If  there  is  nothing  to  compress,  as  in  vacuum,  there 
can  be  no  sound.  The  equations 

^  +  V-(pu)=0  (2.1) 

at 

p[-^r  +  (u-V)u]  +  Vp  =  0  (2.2) 

at 

p  =  p(p)  (2.3) 

characterize  compressible  fluids  where  p  is  the  total  pressure  of  the  fluid,  u  is  a  fluid 
particle  velocity  vector,  and  p  is  the  total  density  of  the  fluid.  Equation  (2.1)  is  a 
mass  conservation  equation,  (2.2)  is  a  lossless  Navier-Stokes  force  equation,  and 
( 2.3)  is  an  equation  of  state.  One  may  add  additional  terms  to  these  equations,  or 
additional  equations  coupled  to  them,  to  account  for  viscous,  heat  conduction,  or 
molecular  relaxation  loss  effects. 

Usually,  in  acoustics,  the  total  variables  p  and  p  do  not  change  very  much 
around  their  ambient  values  p0  and  p0,  and  u  is  assumed  to  have  an  ambient  value 
of  zero,  so  one  may  insert  the  substitutions  p  =  p0  +  p'  and  p  =  p0  +p'  into  the 
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first  two  of  the  above  equations  to  obtain 

<  K  p0  +  p  ) 

- +  V-[(Po  +  //)  u]  =  0  (2.4) 

<n 

(  p0  +  p  ) [  — —  +  ( u-  V )  u]  +  V  ( p0  +  p  )  =  0  ( 2.5) 

ot 

where  //  and  p'  are  the  deviah  ns  of  the  density  and  pressure,  respectively,  from 
their  ambient  values.  Further,  one  can  expand  the  equation  of  state  (2.3)  in  a 
Taylor’s  series11 

P  =P0  +  As  +  |j-s2  +  +  .  .  .  (2.6) 

where  the  condensation  s  =  ( p  —  p0  )/p0  =  p'/p0 . 

The  linear  acoustic  approximation  neglects  any  terms  which  contain  more  that 
one  acoustic  variable  or  its  derivative,  or  power  higher  than  one.  If  one  denotes  the 
variables  u,  p\  and  p'  as  the  acoustic  variables  and  keeps  only  the  firstrorder  linear 
terms,  the  linear  acoustic  equations  result 


where  c0  is  the  speed  of  sound  coming  out  of  the  first  term,  A,  of  the  Taylor’s 
series  expansion  of  the  equation  of  state.  Here,  constant  p0  and  p0  in  time  and 
space  are  assumed.  One  can  solve  these  equations  simply  by  combining  them  into  a 
wave  equation  in  one  of  the  three  variables,  for  example,  the  acoustic  pressure  wave 
equation. 


2 


O' 


2r~  2 


ct 


T-*o* 


P’  -  0 


l 

(2.10) 


Typically,  one  then  Fourier  transforms  away  the  time  dependence  of  these  equations 
to  obtain  the  Helmholtz  equation 


2 


+  k2  p'  —  0 


(2.11) 


where  &  is  the  wavenumber  equal  to  jc/c0 ,  and  u  is  the  rad ’an  frequency. 

This  Helmholtz  equation  is  the  basis  for  most  linear  acoustic  propagation  pro¬ 
grams  such  as  the  PE  and  the  FFP.  The  accuracy  of  these  programs  in  modeling 
the  physics  of  the  fluid  depend,  therefore,  on  how  well  the  linear  acoustic  approxi¬ 
mation  made  in  the  last  paragraph  holds  up.  If  any  of  the  neglected  terms  in  fact 
become  significant  compared  to  the  linear  terms  which  were  kept,  the  linear  acoustic 
model  will  break  down. 


Table  2.1,  for  example,  lists  the  ratio  of  p'  to  p0  for  various  sound  pressure 
levels,  referenced  to  20  /iPa  The  abbreviation  SPL  for  sound  pressure  level  will  be 
used  throughout  this  thesis.  In  Table  2.1  it  is  assumed  that  c0  —  343  m/sec  and 
p0  =  1.21  kg/m3,  values  appropriate  for  air.  It  is  clear  that  for  an  SPL  lower  than 
Table  2.1.  Ratio  of  acoustic  and  ambient  densities  versus  SPL. 


SPL,  dB 

P'/Pn 

60 

1.4x  10"7 

80 

1.4x  10"6 

100 

1.4x  10"6 

120 

1.4x  10~4 

140 

0.0014 

160 

0.014 

180 

0.14 
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120  dB,  i>'  is  negligible  compared  to  />„,  but  for  a  higher  SPL,  especially  above  160 
dB,  p  is  significant  compared  to  p0 . 

Nonlinearity  is  a  cumulative  process.12  A  distortion  of  1/1000  does  not  seem 
appreciable  over  one  wavelength,  but  will  manifest  itself  as  a  serious  distortion  in 
1000  wavelengths. 

To  see  that  nonlinearity  is  cumulative,  look  at  one  of  the  simplest  nonlinear 
acoustic  equations.  The  one-directional  version  of  the  linear  wave  equation  (2.10) 
in  one  dimension  is 


dt  0  dx 


=  o . 


(2.12) 


This  equation  represents  a  plane  wave  moving  in  the  +x  direction.  If  the  second- 
order  nonlinear  terms  but  none  higher  are  included,  the  equation  becomes 


The  plane  wave  now  travels  at  a  speed  c0  1  +  3  ^  '  instead  of  c0 ,  and  this  speed 

W 

is  dependent  on  the  amplitude  of  the  wave.  Therefore,  the  different  parts  of  the 
wave  will  propagate  at  different  speeds  and,  over  long  distances,  the  nonlinear  effect 
can  severely  distort  the  wave  even  though  the  effect  is  small  over  short  distances. 


The  constant  3  is  the  parameter  of  nonlinearity  and  it  contains  two  types  of 
nonlinearity  in  its  equation  3  =  1  +  B/2A .  The  1  comes  from  the  convection  effect 
of  an  acoustic  particle  traveling  at  the  speed  c0  +  |  u  |  instead  of  c0 .  The  B /A  term 
comes  from  the  second-order  term  in  the  nonlinear  equation  of  state  (2.6).  Typical 
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values  of  B/A  are  0.4  and  5.0  for  air  and  water,  respectively,  so  3  is  about  1.2  and 
3.5.  Theories  using  3  account  only  for  plane  wave  nonlinear  effects,  as  there  are 
other  second-order  nonlinear  terms  3  does  not  include.  Nevertheless,  researchers 
have  extended  the  use  of  3  to  the  approximate  cases  of  spherical  waves  and  cylindr¬ 
ical  waves. 

In  the  1960s  Blackstock  developed  a  theory  of  finite  amplitude  sound  that  could 
predict  three  distances  to  describe  the  various  stages  of  weak  nonlinear  propaga¬ 
tion.13,  14  For  an  initially  sinusoidal  wave  Blackstock  called  the  distance  in  which  the 
wave  would  develop  its  first  discontinuity  the  discontinuity  distance.  Later,  when 
the  wave  fully  developed  a  shock,  this  larger  distance  was  called  the  shock  distance. 
Eventually,  the  wave  would  travel  so  far  that  the  attenuation  in  the  medium  would 
begin  to  dominate  the  finite  amplitude  effects,  and  this  was  named  the  old  age  dis¬ 
tance.  For  plane  waves  the  discontinuity  distance,  x,  is 


where  e  is  the  acoustic  Mach  number  |  u  J /c0  and  k  is  the  wavenumber.  For  the 
approximate  case  of  divergent  spherical  waves,  however,  the  discontinuity  distance  r 
is 

F  =  r0eJ7^'  (2-15) 

where  r  is  the  range  from  the  center  of  curvature,  and  r0  is  the  effective  radius  of 
the  spherical  source.  The  plane  wave  result  is  valid  under  the  assumption  that  the 
old  age  distance  xmax  >>  x  where 


Table  2.2.  x  for  peak  pressure  amplitudes  and  frequencies. 


SPL,  dB 

x,  10  Hz 

x ,  100  Hz 

x,  1000  Hz 

120 

32.4  km 

3.24  km 

324  m 

140 

3.24  km 

324  m 

32.4  m 

160 

324  m 

32.4  m 

3.24  m 

180 

32.4  m 

3.24  m 

0.324  m 

Table  2.3.  r  and  rmax  for  peak  pressure  amplitudes  and  frequencies. 
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were  calculated  using  values  of  o  appropriate  for  the  frequencies  involved,  and  the 
values  of  rmax  were  found  using  an  iterative  computer  program.  Muir  and  Carsten- 
sen  have  stated  that  shocks  will  not  develop  if  rmax  <  r  so  Table  2.3  does  not  give 
numerical  values  for  these  cases.  For  each  frequency  and  SPL  xmax  was  always 
greater  than  the  values  of  x  given  in  Table  2.2. 


From  these  tables  it  is  apparent  that  the  discontinuity  distance  decreases  as  the 
frequency  or  the  SPL  increases.  It  is  also  clear  that  the  discontinuity  distances  for 
spherical  waves  are  greater  than  for  plane  waves.  For  example,  at  160  dB  SPL  and 
1000  Hz,  x  is  3.24  m,  but  r  is  25.4  m.  It  makes  sense  that  the  disco ntiniiity  dis¬ 
tance  is  greater  for  the  spherical  case,  since  the  amplitude  of  a  spherical  wave  falls 
off  as  1/r  while  the  amplitude  of  a  plane  wave  changes  only  slightly  when  propagat¬ 
ing.  Unfortunately,  these  results  are  valid  only  for  sinusoids  and  are  approximate 
for  the  cylindrical  and  spherical  wave  cases. 

A  model  that  includes  both  the  parameter  of  nonlkiearity  3  and  absorption  has 
become  known  as  Burgers’  equation.  This  equation  can  model  both  pulses  and 
sinusoids  and  may  be  derived  from  the  form 


dp' 

dt 


+  co 


1  +  3- 


Po  Co 


dp'  _  i  <?y 


dx 


dt‘ 


(2.18: 


where 


H  4_  Ae_  +  (y  -  1)k 

2p0  3  p  Cp 


(2.19) 


accounts  for  the  classical  and  bulk  relaxational  dissipative  processes.  Note  that  the 
only  difference  between  Equations  (2.13)  and  (2.18)  is  the  rightrhand  side.  Here,  p 
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is  the  shear  viscosity  coefficient,  is  a  bulk  viscosity  coefficient,  is  the  ratio  of 
specific  heats,  /•>'  is  the  thermal  conductivity,  and  cp  is  the  coefficient  of  specific  heat 
at  constant  pressure.  Burgers’  equation  is  often  written  as 

Uf  -f-  3uux-  =  flux-,.'  (2.20) 

where  the  plane  wave  particle  velocity  u  —  p'/p0c0  has  been  substituted  into  Equa¬ 
tion  (2.19)  with  the  change  of  variable  x'  =t  —x/c0.  Much  research  has  been 
done  using  Burgers’  equations.15'17 

Trivett  and  Van  Buren  solved  a  Burgers’  type  equation  numerically  via  a  cou¬ 
pled  Fourier  series  analysis  for  plane,  cylindrical,  and  spherical  waves.1819  A  recent 

addition  to  the  study  of  Burgers’  equation  is  given  by  Mitome.20  However,  at  the 
time  of  these  studies  it  was  not  clear  how  molecular  relaxation  effects  could  be 

added  to  the  model,  as  it  is  today.21  Hence,  another  method  of  computation  using 
the  parameter  of  nonlinearity  ,3  arose,  called  the  Pestorius  algorithm. 

The  Pestorius  algorithm  incorporates  the  competitive  effects  of  absorption  and 
finite  amplitude  nonlinearity  by  assuming  the  two  processes  are  independent  over 
small  steps  in  range.22'24  A  small  step  in  range  is  taken  in  the  time  domain  account¬ 
ing  for  nonlinear  steepening.  Then  a  fast  Fourier  transform  is  taken  to  the  fre¬ 
quency  domain  where  losses  are  accounted  for  over  the  range  step.  Then  an  inverse 
fast  Fourier  transform  returns  the  wave  to  the  time  domain,  and  the  process  contin¬ 
ues.  The  iteration  repeats  until  the  modeled  waves  have  propagated  as  far  as  is 
necessary.  Researchers  have  used  this  method  successfully  for  the  prediction  of 
shock  wave  rise  times  from  sonic  booms,  ballistic  waves,  blasts  and  electric  spark 
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pulses.25'27  It  has  been  app'ied  exactly  to  plane  wave  propagation  and  approximately 
to  cylindrical  and  spherical  waves  since  it  employs  the  plane  wave  parameter  of  non¬ 
linearity  3. 

Another  method  based  on  the  parameter  of  nonlinearity  3  is  the  NPE.  which 
stands  for  nonlinear  parabolic  equation.28  This  technique  is  similar  to  the  standard 
parabolic  equation,  but  contains  the  plane  wave  nonlinear  steepening  effects  in  addi¬ 
tion  to  the  effects  of  refraction  and  diffraction.  Here,  one  steps  the  equations  in  the 
time  domain  in  a  particular  direction. 

Aanonsen  et  al.  have  introduced  a  similar  one  way  wave  solution,29'81  this  time 
including  all  the  second-order  nonlinear  effects  coming  from  the  compressible  fluid 
equations,  not  just  those  associated  with  3.  These  authors  solve  this  parabolic  type 
equation  numerically  using  an  implicit  finite  difference  solution  of  a  frequency  cou¬ 
pled  Fourier  series  expansion,  not  unlike  the  method  of  Trivett  and  Van  Buren. 

The  Aanonsen  et  al.  work  seems  most  promising,  but  like  the  NPE,  the  Pes- 
torius  algorithm,  and  Burgers’  equatiun  methods,  it  can  model  only  waves  which  are 
propagating  primarily  in  one  direction.  These  methods  are  not  sufficient  to  model 
the  propagation  of  sound  outdoors,  near  ground  surfaces  which  reflect  sound. 
Reflected  waves  are  not  accounted  for  in  the  above  theories. 

In  addition,  only  the  Aanonsen  method  accounts  for  all  the  second-order  non- 
linearities  of  the  compressible  fluid  equations.  The  second-order  nonlinearities  not 
associated  with  3  are  important  in  the  near  field  of  a  blast  source,  for  example, 
since  the  waves  bene  are  far  from  planar.  Second-order  terms  are  also  important 
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when  direct  and  reflected  waves  interact,  near  a  surface  or  near  a  turning  point  in  a 
refractive  atmosphere.  Chapter  4  will  give  a  model  to  include  all  the  second-order 
nonlinear  effects  and  both  incoming  and  outgoing  waves  simultaneously. 


3.  SOUND  PROPAGATION  OUTDOORS 


This  chapter  reviews  the  keys  to  the  modeling  of  sound  propagation  outdoors. 
Much  of  this  matenai  has  appeared  as  summary  articles  in  the  literature'52"34  and  is 
widely  accepted.  Other  than  geometrical  spreading,  the  primary  mechanisms  which 
can  influence  tne  prediction  of  sound  levels  over  flat,  featureless  surfaces  are  ab¬ 
sorption,  ground  impedance,  refraction,  and  turbulence. 

The  best  understood  aspect  of  sound  propagation  in  the  free  held  is  the  absorp¬ 
tion  and  transference  of  energy  as  an  acoustic  wave  moves.  The  primary  processes 
are  viscous  and  heat  conduction  effects,  and  molecular  relaxalional  effects.  The 
classical  effects  were  discovered  in  the  1800s,  whereas  the  understanding  of  relaxa¬ 
tion  has  become  clear  only  in  the  last  25  years.  This  thesis  will  not  detail  either  of 
these  processes,  as  a  number  of  standard  textbooks  describe  the  particulars  of  the 
mechanisms. 36  • 36 

The  theory  to  account  for  the  influence  of  an  outdoor  ground  surface  is  well 
understood  also,  if  the  impedance  Z  of  the  ground  is  known.  Two  of  the  better 
journal  articles  which  review  the  important  aspects  of  the  development  of  this 
theory  are  listed  in  the  references,37'3**  although  the  latter  has  innumerable  mis¬ 
prints.  Unfortunately  this  has  been  a  problem  in  the  literature  of  propagation  near 
finite  acoustical  impedance  boundaries,39  the  difficulty  often  being  that  some  au¬ 
thors  use  the  notation  and  others  e  ~‘~t.  In  addition,  Nobile  and  Hayek  recent¬ 
ly  have  found  an  exact  solution  for  sound  propagating  in  a  homogeneous  atmo¬ 
sphere  with  arbitrary  impedance  Z  that  is  easily  implemented  in  a  computer  pro- 
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gram.40  Such  a  program  will  be  used  in  Chapter  11  to  examine  what  peak  levels 
result  in  the  linear  acoustic  limit  for  certain  frequencies  and  impedances,  for  com¬ 
parison  with  the  nonlinear  predictions  described  there. 

The  problem,  however,  is  that  researchers  have  had  difficulty  accurately  ascrib¬ 
ing  surface  impedances  to  particular  ground  surfaces.  Numerous  models  have  been 
published  in  the  literature  to  relate  the  physical  properties  of  the  ground  to  Z  as  a 
function  of  frequency.  One  of  the  more  frequently  cited  theories  is  provided  by 
Chessel,  based  on  the  work  of  Delany  and  Bazley.41  In  this  formulation  the  ground 
is  characterized  by  a  static  flow  resistance,  a  measure  of  how  easily  air  moves 
through  a  ground  sample.  The  Delany- Bazley-Chessel  theory  is  simple  and  is  com¬ 
monly  in  use  in  most  modem  outdoor  acoustic  propagation  programs  such  as  the 
FFP  (fast  field  program)  and  the  PE  (parabolic  equation  method).  Another  more 
complex  model  of  natural  outdoor  surfaces  is  the  four  parameter  model  of  Atten¬ 
borough.  This  theory  includes  a  second  physically  intuitive  parameter  in  addition  to 
the  flow  resistance,  called  the  porosity.  The  porosity  is  the  volume  fraction  of  air  in 
a  ground  sample.  In  addition  to  porosity  and  flow  resistance,  Attenborough  adds 
two  more  parameters,  the  tortuosity  and  the  pore  shape  factor  ratio.42  The  tortuosity 
is  a  measure  of  how  tortuous  or  twisty  the  pores  of  the  ground  are,  and  the  shape 
factor  ratio  gives  some  measure  of  the  pores’  cross-sectional  shape  relative  to  a 
square  ora  circle. 

It  will  become  clear  in  Chapter  4  that  the  second-order  nonlinear  equation  of 
the  air  can  be  formulated  only  in  the  time  domain.  All  of  the  above  cited  ground 
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models,  on  the  other  hand,  are  in  terms  of  the  impedance  Z  in  the  frequency 
domain.  Thus,  an  efficient  solution  to  the  outdoor  propagation  problem  will  require 
a  formulation  of  the  ground  in  the  time  domain.  Morse  and  Ingard  have  described 
such  a  model  for  propagation  in  a  porous  medium.43  In  addition,  Attenborough  has 
related  this  theory  to  a  low  frequency  model  of  the  ground.  Chapter  5  will  elaborate 
on  this  model  m  detail. 

Morse  and  Ingard’ s  equations  do  not  include  the  effects  of  nonlinearities  in  a 
porous  ground,  however.  Little  work  has  been  performed  on  this  topic.  One  excep¬ 
tion  has  been  the  research  by  Kuntz,  Nelson,  and  Blackstock.44  46  These  gentlemen 
investigated  the  effects  of  nonlinear  propagation  in  batted  Kevlar,  a  strong  fibrous 
material.  Their  results  showed  that  the  nonlinearity  in  porous  propagation  is  associ¬ 
ated  with  the  flow  resistivity  and  is  not  similar  to  the  nonlinearity  in  air.  Hersh  also 
developed  a  model  for  the  acoustic  properties  of  materials  such  as  Kevlar.46  Wilson, 
McIntosh,  and  Lambert  have  extended  some  of  the  knowledge  in  this  area,47  but 
more  research  is  needed.  Note  that  the  porous  model  employed  in  Chapter  5  is  a 
linear  model. 

Another  m^jor  factor  in  predicting  the  propagation  of  sound  outdoors  is  refrac¬ 
tion  in  the  atmosphere.  Because  of  the  heating  and  cooling  of  the  ground  and  be¬ 
cause  of  the  effect  of  wind,  the  speed  of  sound  in  the  air  varies  as  a  function  of 
height,  distance,  direction,  and  time.  Simple  refraction  models  have  been  incor¬ 
porated  into  programs  such  as  the  FFP  and  the  PE  to  estimate  the  combined  effects 
of  ground  impedance  and  refraction. 
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Random  inhomogeneities  in  the  air,  called  turbulence,  also  affect  the  propaga¬ 
tion  of  sound  outdoors.48,49  This  phenomenon  is  particularly  important  for  long 
ranges,  because  the  weak  scattering  of  acoustical  energy  by  turbulence  can  accumu¬ 
late  over  distance.  The  uncertainty  of  turbulence  effects  renders  all  acoustic  propa¬ 
gation  prediction  programs  approximate,  and  more  research  is  needed  on  this  topic. 
Some  recent  simulation  results  and  references  to  the  literature  may  be  found  in 
McBride.60 

This  dissertation  later  develops  a  model  which  at  present  accounts  for  neither 
refraction  nor  turbulence.  Refraction  should  be  no  problem  to  add  at  a  later  date, 
as  long  as  one  rederives  the  equations  in  the  next  chapter  under  the  assumption  that 
p0  and  c0 ,  the  ambient  density  and  speed  of  sound,  are  functions  of  space.  There 
should  be  no  need  to  regard  them  as  functions  of  time  since  they  change  negligibly 
during  the  passage  of  a  sound  wave.  One  might  implement  turbulence  as  a  pertur¬ 
bation  on  the  variables  p0  and  c0  as  functions  of  space,  making  many  runs  with 
many  different  perturbations  to  give  some  statistical  average  propagation  results. 
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4  THE  TIME  DOMAIN  AIR  PROPAGATION  EQUATIONS 

This  chapter  presents  the  equations  to  describe  lossy,  weakly  nonlinear  sound 
propagation  in  air.  The  equations  include  both  of  the  classical  dissipation  effects, 
the  effect  of  a  bulk  viscosity,  and  all  of  the  second-order  nonlinear  terms  from  the 
compressible  fluid  equations.  The  expressions  do  not  include  relaxation  effects. 
The  equations  are51'63 


The  variables  introduced  in  Chapter  2  are  used  in  these  equations,  except  for  p, 
which  now  represents  the  acoustic  pressure  deviation  instead  of  the  total  pressure. 
In  Chapter  2  p’  denoted  the  acoustic  pressure.  Again,  these  equations  assume  p0 
and  pQ  are  constants  in  time  and  space. 

The  newly  introduced  variables  are  the  following:  T\  the  acoustic  temperature 
deviation;  Sfr,  the  acoustic  entropy  deviation  with  vibrational  degrees  of  freedom 
frozen;  T0,  the  ambient  temperature;  B/A,  the  ratio  of  the  first  two  terms  of  the 
Taylor's  series  expansion  of  total  pressure  in  terms  of  p';  and  0,  the  coefficient  of 
thermal  expansion.  Note  that  in  the  rest  of  this  thesis  0  does  not  represent  the 
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parameter  of  nonlinearity  as  it  did  in  Chapter  2.  For  ideal  gases,  3T0  =  1.  The 
Equations  (4.1)  to  (4.5)  use  the  notation  of  Pierce,  with  the  second-order  nonlinear 
terms  of  Hamilton  and  Blackstock. 

To  clarify  the  above  equations,  one  can  write  them  as 
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(4.15) 

(4.16) 

(4.17) 


Equations  <  4.6)  to  (4.10)  give  the  linear  lossless  terms  explicitly.  The  nonlinear 
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terms  begin  with  the  letter  N  in  (4.11)  -  (4.13),  and  the  dissipative  terms  begin  with 
the  letter  D  in  (4.14)  -  (4.17).  If  N 1,  N 2,  and  N3  — *-0,  Equations  (4.6)  -  (4.10) 
become  linear  and,  similarly,  if  Dl,  D  2,  03,  and  04  — *-0,  (4.6)  -  (4.10)  become 
lossless. 

By  examining  the  terms  Nl,  N 2,  and  JV 3,  one  may  see  that  each  of  these 
terms  is  indeed  second-order.  They  are  each  made  of  the  sums  or  differences  of 
two  acoustic  quantities  or  their  derivatives  multiplying  each  other.  In  the  case  of 
N 3  the  acoustic  density  deviation  is  multiplied  by  itself  to  give  p’2.  Note  that  there 
are  no  higher-order  quantities  such  as  the  product  of  three  acoustic  variables  or 
their  derivatives. 

Since  these  terms  contain  products  of  acoustic  variables,  it  is  impossible  to 
Fourier  transform  Equations  (4.6)  to  (4.10)  without  introducing  convolution,  a 
computationally  expensive  process.  This  circumstance  is  not  unique  to  the  lossy 
nonlinear  acoustic  equations  employed  here,  but  is  a  mathematical  reality  for  any  set 
of  nonlinear  partial  differential  equations.  Hence  this  research  has  centered  on  time 
domain  solutions. 

To  make  a  time  domain  solution  possible  it  is  necessary  that,  each  of  the  Equa¬ 
tions  (4.6)  -  (4.10)  have  no  more  than  one  time  derivative  present  The  only 

offender  is  (4.7)  which  has  two.  Here,  the  time  derivative  term  p0— —  must  be 

at 

present  regardless  of  the  amplitude  of  the  wave.  But  since  the  term  p'— —  in  N2  is 

at 

not  present  in  the  linear  limit  one  can  eliminate  this  time  derivative  as  follows: 
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N2  =  -p'^-p0wVu 

CJt 

Po 

=  Vp  -  p0u-Vu. 
Po 


-  p0u-V u 


(4.18) 


Here,  the  linear  lossless  portion  of  ( 4.7T  has  been  substituted  into  the  term  N 2  it¬ 
self.  Note  that  the  new  form  of  N 2  is  still  of  second-order. 

Before  a  particular  geometry  is  specified,  the  following  assumption  should  be 
understood.  Equations  (4.1)  -  (4.5),  or  equivalently  (4.6)  -  (4.18),  implicitly  as¬ 
sume  that  the  interaction  effects  of  dissipation  and  second-order  nonlinearity  are 
small.  If  one  were  to  include  such  effects,  terms  would  be  introduced  that  were 
third-order  or  higher.  One  neglects  such  interaction  effects  to  keep  the  equations 
consistently  second-order. 

Note  that  it  is  easy  to  incorporate  the  effects  of  molecular  relaxation  into  the 
above  equations.  To  include  the  primary  relaxation  processes  of  oxygen  and  nitro¬ 
gen  in  air.  Equation  (4.3)  becomes 

dsfr  dT  i  dTo  0 

PoT°  +  p°Cv'~dT  +  PoC^~dT  -  kV  r  (419) 

where 


and 


dTi  T  _  Ti 
dt 


(4.20) 


ar  2  t  -  t2 

dt  T  2 


(4.21) 
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The  subscripts  1  and  2  denote  the  relaxation  processes  of  oxygen  and  nitrogen, 
respectively.  Using  the  notation  of  Pierce,  Tx  and  T2  are  apparent  vibration  tem¬ 
peratures,  c0i  and  Cy2  are  specific  heats  at  constant  volume  associated  with  internal 
molecular  vibrations,  and  rx  and  r2  are  relaxation  times  for  vibrational  energies. 
The  introduction  of  Equations  (4.19),  (4.20),  and  (4.21)  poses  no  additional 
difficulties.  One  merely  needs  to  keep  track  of  two  additional  variables,  Tx  and  T2. 

Now  specialize  the  Equations  (4.6)  -  (4.18)  for  a  two-dimensional  cylindrical 
geometry,  natural  for  outdoor  sound  propagation.  This  thesis  will  use  the  coordi¬ 
nates  ( d ,  z,  <d),  as  seen  in  Figure  4.1.  Here,  the  variable  d  represents  range,  in¬ 
stead  of  r  because  this  is  a  spherical  coordinate  variable,  and  instead  of  p  because 
one  might  confuse  this  with  a  density  variable. 

Tf  there  are  no  variations  with  iv.opcct  to  the  variable  <f>,  i.e.  — —  =  0  and  the  6 

d<p 

A 

component  of  u  is  zero,  then  u  =  ud  +  ve  and  Equations  (4.6),  (4.7),  and  (4.8)  be¬ 
come 


/  \ 


dp'  du  ,  u  ,  dv  ATi  a 

dt  +P°  dd  +  d  +  dz  =NU 

(4.22) 

dry  ■»  ▼  * 

dt  +  dd  24  +  DL4 

(4.23) 

P0~  +  ^  =  N2B  +D1B 
dt  dz 

(4.24) 

dSfr 

P„T0  £  =D2A 

(4.25) 

where 


and 


NIA  =  -p' 


du  u  dv 
dd  d  dz 


dp '  dp 

-  u-h —  v- 


dd 


dz 


N2A  =  -  p0 

Po  dd 


8u  ,  du 

U - h  V - 

dd  dz 


N2B  =  -  p0 

Po  dz 


dv  ,  dv 

U - h  V - 

dd  dz 


D  1A  —  (  —fl  +  Pg  ) 
o 


cfiu  1  du 


u  cPv 


+  p 


d?u 


+ 


dd2  d  dd  d2  dddz 

1  du  d2u  u 

+  — 7T - itt 


dd2  d  dd  dz2  d2 


DIB  =  ( — /i  +  p.B) 


(Pu  1  du  d?v 

dddz  d  dz  dz 2 


+  P 


cPv  ,  1  dv  ,  d*v 

•  1  r\  1  I 


dd2  ddd  dz2 


D2A  s*  k 


dPT'  ,  1  dT'  &T 


dd2  d  dd 


dz 2 
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(4.26) 

(4.27) 

(4.28) 


(4.29) 


(4.30) 


(4.31) 


Equation  (4.22)  simplifies  under  the  assumption  of  a  homogeneous  atmosphere. 
For  this  case 


thus 


dp0  =  dpo_ 
dd  dz 


(4.32) 


dp 

dt 


-  +  is  [(^  +  "'H  +  la  [h  + '’’H  —  [p° +  P'H  ■  (433) 


Further,  one  can  make  rearrangements  and  see  that  the  four  time-dependent  equa¬ 


tions  can  take  the  forms 
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dp_ 

dt 


—  +  4-  —  =  N2D  +  DID 

dt  dz  p0 


where 


dSfr 

dt 


=  D2B 


(4.34) 

(4.35) 

(4.36) 

(4.37) 


NIB  =-(p0  +  p')^  (4.38) 

and  N2C  =N2A/p0,  N2D  =  N2B/p0,  DIC  =  DlA/p0,  DID  =  D\B/p0,  and 
D2B  =  D2A /p0  T0 . 


Equations  (4.34)  to  (4.37)  have  been  written  in  a  special  way.  They  are  each 
almost  written  in  conservation  law  form  where  the  four  acoustic  variables  p\  u,  v, 
and  Sfr  are  the  only  time  derivatives  which  occur.  The  terms  NIB ,  N2C ,  N2D, 
DIC,  DID,  and  D2B  may  be  regarded  as  “pseudo  -  source”  terms.  Also  note 
that  when  these  right-hand  side  terms  approach  zero  and  p'  becomes  negligible  com¬ 
pared  to  p0 ,  the  equations  accurately  describe  two-dimensional  linear  acoustics.  In 
contrast,  if  one  uses  the  full  equations  of  computational  fluid  dynamics  to  solve  a 
nearly  linear  acoustics  problem,  the  numerical  solutions  are  usually  poor  in  the 
linear  limit  Now  one  can  write  these  lossy,  second-order  nonlinear  acoustical  equa¬ 
tions  as  the  system 


wt  +  Fj  -t-  Gz  =  SF  ■+■  SG 


where 


( 4.39) 
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and 


(4.40) 


(4.41) 


(4.42) 


(4.43) 


0 

***  ~  N2D  +  DID 
0 


(4.44) 


with  the  equations  of  state  (4.9)  and  (4.10).  The  system  (4.39)  is  a  form  solvable 
numerically  by  time  stepping.  Given  the  initial  conditions  of  p\  u,  v,  and  Sfr  over 
all  space,  these  variables  and  p  and  T'  are  known  everywhere  at  all  future  times  us¬ 
ing  (4.39)  and  the  equations  of  state  (4.9)  and  (4.10).  Chapter  7  will  present  a 
specific  numerical  solution  method  for  (4.39). 
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5.  THE  TIME  DOMAIN  POROUS  MEDIA  EQUATIONS 

This  chapter  examines  the  time  domain  phenomenological  equations  of  sound 
propagation  in  a  porous  medium.  The  porous  equations  are  compatible  with  the 
time-dependent  equations  of  the  air.  The  porous  medium  equations  will  be 
analyzed  in  the  frequency  domain  to  determine  the  surface  impedance,  effective 
sound  speed,  and  absorption  of  a  porous  ground  surface  model  based  on  the  equa¬ 
tions.  Chapter  9  will  describe  the  time  domain  numerical  solution  used  for  the 
porous  medium. 

Sound  propagation  in  pores  is  a  classical  topic,  which  has  been  examined  since 
the  beginning  of  research  in  acoustics.  In  1896  Lord  Ravleigh  devoted  several  sec¬ 
tions  in  his  book,  The  Theory  of  Sound,  to  the  propagation  of  sound  in  pores  and 
narrow  tubes.54  Rayleigh’s  work  was  based  on  research  by  Kirchhofif.  The  book  that 
is  probably  most  often  cited  on  porous  materials,  however,  is  Sound  Absorbing  Ma¬ 
terials  by  Zwikker  and  Kosten  from  1949.66  They  extended  Rayleigh’s  calculations  in 
pores  to  acoustical  absorbing  materials.  Zwikker  and  Kosten  include  many  expen- 
mental  results  along  with  their  theoretical  analysis. 

This  chapter  employs  a  model  that  is  originally  from  Morse  and  Ingard’s 
Theoretical  Acoustics .56  These  authors  model  a  bulk  porous  medium  as  an  idealized 
fluid.  Their  model  porous  acoustical  equations  are 

n^-  +  p0V'U  =  o  (5.D 


Pp  ~  +  <*>«  +  \Tp'  =  0 


(5.2) 
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P  =  KpPoP’ 


(5.3) 


where  f!  is  the  porosity,  p'  and  pn  are  the  acoustic  and  ambient  densities,  u  is  the 
average  particle  velocity  in  the  direction  of  propagation,  pp  is  an  effective  density,  4> 
is  a  flow  resistance,  p'  is  the  acoustic  pressure,  and  Kp  is  the  effective  compressibili¬ 
ty  of  the  fluid  in  the  pores.  Equations  (5.1),  (5.2),  and  (5.3)  are  the  equations  of 
continuity,  motion,  and  state,  respectively.  Note  that  k  and  pn  should  be  viewed 

r  r 

as  functions  of  frequency.  In  fact,  Morse  and  Ingard  solve  (5.1)  -  (5.3)  in  the  fre¬ 
quency  domain,  bypassing  the  time  dependence  of  these  equations.  Their  transfor¬ 
mations  using  e~l~(  time  dependence  were 


Pe  =  Pp  1  +  - 

( 5.4) 

Pp M 

1 

c'  y/o.*. 

( 5.5) 

C 

II 

( 5.6) 

for  an  effective  density,  sound  speed,  and  compressibility  to  obtain  the  Helmholtz 
equation 


( 5.7) 


where  u  is  calculated  from  a  =  —V  xl\  The  impedance  of  such  a  porous  material  is 
then 


Z  -  pece  =  (pe/Ke)% 

=  (Pp/K-pft  )*(  1  +  i $/pp *j) *  . 


(5.8) 


Using  Morse  and  Ingard’ s  results,  but  using  the  ewt  time  dependence,  Donato  com- 
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pared  this  impedance  modei  to  some  averaged  experimental  data  over  a  wide  range 

of  frequencies  in  1977. 57  He  showed  that  the  above  model  was  probably  not  as  good 
as  a  model  where  porosity  varied  as  a  function  of  depth. 

Attenborough  has  analyzed  porous  equations  similar  to  those  of  Morse  and  In- 
gard  and  then  related  them  to  a  low  frequency  model  of  a  finite  ground 
impedance.58'60  His  form  of  the  generalized  pore  equations  in  one  dimension  are 


ni2-  +  ,4“  =o 

dt  p  dx 


K  du 


n  dt  +<f,“ +  If  =  0 


Po 

C02(p  ~  Po')  ^  1<P  - Po > 


(5.9) 

(5.10) 

(5.11) 


Note  that  (5.9)  was  misprinted  in  Attenborough’s  1983  J.  Acoust.  Soc.  Am.  paper, 
but  appeared  correctly  in  his  1982  Physics  Reports  work.  These  equations  use  a 
cmTerent  notation  than  (5.1),  (5.2),  and  (5.3).  Here  p  represents  a  total  pressure, 
with  p0  the  ambient  pressure,  and  similarly  for  p  and  p0.  The  constant  K 
represents  an  effective  density  factor  or  structure  factor.  If  one  linearizes  these 
equations  to  first-order,  under  the  assumption  p  —  p0  +  p'  and  p  =  p0  +  p\  and 
remembers  that  at  low  frequencies61  -y  — ►  1, 


fi 


Po 


dp’  ,  n  dn  a 
at  +p°ax  =0 

(5.12) 

du  ,  dp'  n 

—  +  Tm  +  -f-  =  0 
at  dx 

(5.13) 

II 

(5.14) 
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These  equations  involve  only  linear  variables  and  may  be  compared  to  Morse  and 
Ingard’s  Equations  (5.1)  -  (5.3). 

To  derive  equations  for  a  complex  wavenumber,  one  eliminates  p'  and  p’  to  ob¬ 
tain  the  governing  wave  equation  for  this  porous  model  in  the  average  velocity  u, 


v  ,  a.  n  du  2  d2u  A 


(5.15) 


By  substituting  u  =  e 
results, 


i(*jt  -  kxx )  ■ 


into  this  equation,  the  plane  wave  dispersion  relation 


or 


—<jj2p0K  +  n  +  k2p0c 2  =  0 


(5.16) 


.*n_ 

Po^o 


■U)  . 


(5.17) 


This  latter  form  will  be  useful  later,  since  the  effective  sound  speed  in  the  pores  is 
ceff  =  nj/Re  [kx  ]  and  the  attenuation  is  a  —  —Im[kx  ],  i.e.,  kx  =  u ;/ceff  —  ia. 

In  his  1983  paper  Attenborough  relates  the  structure  factor  K  to  one  of  the 
parameters  of  his  four-parameter  model  of  a  porous  medium  by 


K 


(5.18) 


where  q  is  the  tortuosity.  The  correspondence  between  the  Morse  and  Ingard 
model  and  his  four-parameter  model  is  exact  when 

n2 

—  =  1  (5.19) 

s 


where  n  is  the  dynamic  shape  factor  and  s  is  the  static  shape  factor.  In  the  four- 
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parameter  theory  of  this  1983  paper,  n  =  s  =  1  corresponds  exactly  to  the  case  of 
perfectly  circular  cylindrical  pores.  Note  that  in  his  1987  paper  Attenborough  shifts 
to  the  notation  sp  =  n2A2/s  where  A  =  1/2  for  circular  cylindrical  pores,  where  sp 
is  called  a  pore  shape  factor  ratio.  Under  the  1983  theory,62  a  typical  value  of  q  for 
soil  is  1.5. 


In  the  1983  paper  the  impedance  of  the  Morse  and  Ingard  phenomenological 
equations  is  found  to  be 

— —  =  (2~,f2  )~%n(  1  +  i)(^/p0sui)%  (520) 

Poco 


which  is  independent  of  tortuosity.  Note  that  the  exponent  of  (2yfl )  is  misprinted 
in  the  1983  paper  itself,  but  (5.20)  may  be  verified  from  his  Equations  (32),  (36), 
and  ( 38)  under  the  assumption  of  large  $  and  low  c o. 

For  a  typical  outdoor  mown  grass  surface  $  and  fi  might  be  300,000  mks 
Rayls  =  300,000  kg/(s  m2)  and  0.30,  respectively.63  For  these  values  and  where 
p0  =  1.21  kg/m3,  ->  — *•  1  (low  frequency),  and  s  =  n  —  1  (circular  cylindrical  pores) 


Z 

Poco 


256.44 


( 1  +  i )  . 


(5.21) 


This  is  the  impedance  predicted  from  the  Attenborough  interpretation  of  the 
phenomenological  porous  model.  Table  5.1  lists  values  of  Z/tp0c0(l  +  t ) 3  for  vari¬ 
ous  frequencies,  using  the  above  constants.  With  these  same  values,  one  can  find 
the  speed  of  sound  and  attenuation  in  the  pores  from  (5.17).  Here, 
kx2  =  25.5 t,imes  10-6^2  -  t0.632u\  Column  three  in  Table  5.1  gives  values  of  the 
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Table  5.1.  Porous  model  impedance,  sound  speed,  and  absorption. 


f,  Hz 

Z/p0c0(l+i) 

ceff  =  jj/Re[kx],  m/s 

a  =  —Im[kx],  Np/m 

500 

11.468 

93.61 

29.58 

400 

12.82 

84.79 

26.79 

300 

14.80 

74.35 

23.5 

200 

18.13 

61.48 

19.43 

100 

25.64 

44.03 

13.91 

75 

29.61 

38.25 

12.08 

50 

36.26 

31.35 

9.90 

25 

51.28 

22.25 

7.02 

10 

81.09 

14.08 

4.45 

5 

114.68 

9.97 

3.14 

effective  sound  speed  in  the  pores,  and  column  four  gives  the  attenuation  coefficient 
for  the  listed  frequencies. 


Notice  that  as  the  frequency  decreases,  the  magnitude  of  impedance  increases, 
the  effective  sound  speed  decreases,  and  the  attenuation  decreases.  Also  note  that 
below  50  Hz  ceff  is  less  than  l/10th  the  speed  of  sound  in  air,  c0  =  343  m/s.  This 
fact  supports  the  hypothesis  that  the  porous  medium  is  locally  reacting,  i.e.,  plane 
waves  in  air  incident  on  the  porous  surface  at  significant  angles  of  incidence  from 
normal  will  propagate  into  the  pores  normal  to  the  surface.  If  the  surface  were  not 
locally  reacting,  waves  incident  on  one  part  of  the  surface  would  affect  the  waves 
propagating  on  nearby  parts  of  the  surface.  Local  reaction  implies  that  each  part  of 
a  porous  surface  is  independent  of  all  others.  This  concept  is  useful  and  will  be 
used  in  Chapter  9. 

Because  the  reactance  of  the  impedance,  the  imaginary  part  of  Z,  is  equal  to 
the  resistance,  the  real  part  of  Z  in  Equation  (5.20),  surface  waves  should  not 
develop  for  propagation  over  this  porous  medium.  Surface  waves  are  independent 
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propagating  waves  at  the  interface  of  the  air  and  porous  medium  which  can  some¬ 
times  develop  depending  on  whether  the  reactance  of  the  impedance  is  greater  than 
the  resistance.  Surface  waves  have  been  studied  in  great  detail,  and  are  reasonably 
well  understood.64"68  For  a  spherical  source  they  propagate  along  the  ground  with  a 
r~ %  dependence  and  decrease  exponentially  with  height.  Because  these  waves  de¬ 
crease  with  distance  as  r~^,  their  contribution  to  the  total  sound  field  can  be 
significant  at  medium  ranges  since  the  rest  of  sound  field  components  fall  off  as  1/r. 
Surface  waves  are  different  from  the  “ground  waves”  of  AM  radio  propagation,  and 
one  may  obtain  a  more  complete  discussion  on  surface  waves  in  the  references. 
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6.  MODEL  INITIAL  CONDITIONS  IN  AIR 

This  chapter  presents  the  initial  conditions  for  both  blast  waves  and  spark 
pulses  used  in  the  numerical  simulations  given  in  Chapters  10  and  11.  Blasts  are 
discussed  first,  and  then  sparks.  Scaling  law  theory  is  also  discussed  with  reference 
to  the  blast  initial  conditions.  This  chapter  closes  describing  how  the  initial  condi¬ 
tions  were  implemented  near  surfaces. 


6.1  Blast  Initial  Conditions  and  Scaling  Laws 


The  pulse  shape  originally  used  in  this  investigation  was  the  ideal  blast  pulse  of 
Reed.69  The  shape  closely  approximates  the  pressure  waveforms  emitted  both  by 
electric  sparks  and  by  explosions,  on  different  scales.  Reed’s  pulse  is 


Pit)  =  Ap0 


(6.1) 


for  0  <  t  <  t  and  pit)  =0  otherwise.  The  constant  £>p0  is  the  peak  pressure  of 
the  pulse,  r  is  the  waveform  duration,  and  t+  is  the  length  of  the  positive  phase, 
when  pit)  >  0. 

One  problem  with  this  idealization  for  numerical  computations  is  that  it  has  a 
discontinuity  at  t  =  0,  where  it  jumps  from  piO~)  =  0  to  p(0+)  =  Ap0 .  To  elim¬ 
inate  this  discontinuity,  which  contains  much  high  frequency  energy,  a  modified 
Reed  pulse  was  developed.  It  has  the  wave  shape  of  (6.1),  with  a  raised  cosine  in¬ 
troduction  so  that  at  no  point  in  the  pulse  is  pit)  or  p' it)  discontinuous.  The  ex¬ 
pression  for  this  pulse  is 


where 


P\(t),  0  <  t  <  1 3 

pit)  =  p2(t) ,  £3  <  t  <  t2  +  r 

0,  otherwise 


Ap0  L  it  t 

px{t)  =  ~Y~  1  -  cos— 


(6.2) 


t  t  o  t  —  to  t  —  to 

p2(t)=lp0  1-—-L  1 - 1 - -i. 


(6.4) 


Here  tl<  t2<  1 3  where  t2  is  the  “delay”  of  the  original  Reed  pulse  and  where  t: 
is  when  the  two  waveforms  both  meet 


P\(t)  =  p2(t)  t_ 


(6.5) 


and  their  derivatives  match 


4tPi(*)  =  4rP2^) 

t  -  h 


Because  this  thesis  considers  only  weak  nonlinear  acoustical  waves,  the  pulse 
shape  should  be  realizable  in  the  linear  acoustics  limit.  In  this  linear  limit  the  pulse 
pit)  additionally  should  balance,  satisfying  the  conservation  of  impulse 


1 2+  r 

f  pit)  dt  =  0 
o 


(6.7) 


This  expression  means  that  after  the  pulse  passes  by,  it  should  leave  the  compressi¬ 
ble  fluid  particles  in  their  original  state. 
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For  the  blast  wave  case  the  durations  r  =  0.04  s  for  the  Reed  waveform  and 
t  j  =  0.5x  10-3  s  for  the  first  n  radians  of  the  cosine  introduction  were  used.  These 
values  mean  the  0%  to  100%  rise  time  of  the  cosine  intm  diction  was  approximately 
l/80th  the  length  of  the  rest  of  the  pulse.  This  duration  of  r  corresponds  to  that 
observed  30  m  from  a  0.57  kg  C-4  plastic  explosive  detonation.70  The  peak  level  for 
this  pulse  is  164.1  dB  SPL  at  30  m. 

A  computer  program  was  written  to  find  t+,  t2 ,  and  t3  to  satisfy  (6.5),  (6.6), 
and  (6.7)  with  the  above  given  values  of  tx  and  t.  The  values  found  were  t+  = 
11.0326x  10-3,  t2  =  0.50293016x  10~3,  and  t3  =  0.50585915x  10~3  s. 

The  above  specification,  however,  is  only  on  the  pressure.  The  system  of  equa¬ 
tions  (4.39)  requires  initial  values  on  the  four  variables  p\  u,  u,  and  Sfr,  which  im¬ 
ply  p  and  T'.  In  this  research  the  approximation  was  made  that  the  initial  condi¬ 
tions  would  be  lossless  and  linear,  under  the  assumption  that  the  variable  values 
would  settle  into  their  lossy  nonlinear  relationships  in  a  few  time  steps. 

Using  this  lossless  linear  assumption,  Sfr  was  assumed  zero  everywhere  initially. 
Then  T'  is  linearly  related  to  p  from  Equation  (4.5).  Also,  in  a  linear  approxima¬ 
tion  one  can  find  p  =  c02p',  and  u  and  v  from  p  using  a  potential. 


For  spherical  pulses  the  radial  component  of  particle  velocity  vr  is  related  to  the 
pressure  by  the  spherical  potential  Fit  —  r/c0)  as 


at 


q  Fit  -  r/c0) 


(6.8) 


and 


r 
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p  Fit  —  r/c0 ) 

vr - 1 - r} - 

Poco  P0rl 


The  variables  u  and  v  are,  of  course,  related  to  vr  by  the  geometrical  expressions71 
u  =  ursin#  and  v  —  vrcos 8.  Thus,  given  p  for  a  spherical  pulse,  p',  u,  v,  Sfr,  and 
T'  can  be  found  in  a  linear  lossless  approximation. 

Corresponding  to  the  pressure  waveform  (6.2)  -  (6.4)  one  can  find 

Fit  —  r/c0  )/r  by  integrating p  with  respect  to  t.  The  result  is 

Fi  it  —  r/c0  )  =  tb ) 


^>o  ,  h  . 

— —  th - sin  - 

2  b  7T  \tl 


t  q  —  — sin 


0  5:  tb  —  ^3 


+  ^1*  *3  <  tb  ^  t2  +  T 
otherwise 


(6.10) 


where 


[l  I  —  ?2)  —(t3  —  t2)‘ 

h-t3-  7+  —  - 5 - 


1  1  <%  -«2>3-»3-«2>: 


'+  T 


1 1  ,  i  ytb  - *2)4  -  it3  -  t2r 

+  r3  +  t  T2  4 


1  (t[j  —  f2)6  —  (t3  —  £2){ 


(6.11) 


For  all  the  runs  of  Chapters  10  and  11  an  approximate  version  of  Ix  was  used. 
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approximate  I  x 
Mj 


1  ( tb-h >6 
tH+  5 


t^T‘ 


Ub 


t q)4 


-(6.12) 


For  the  particular  values  of  flt  t2>  *3»  f+>  r  used  in  this  study,  the  difference 
between  (6.11)  and  (6.12)  is  very  small.  Figure  6.1  plots  Equation  (6.10)  using 
both  Equations  (6.11)  and  (6.12),  with  the  former  using  a  solid  line  and  the  latter  a 
dashed  line. 

A  most  successful  theory72  for  predicting  the  effects  of  explosions  has  been  the 
scaling  laws  of  Hopkinson  from  1915.  These  relationships  have  been  experimentally 
verified  over  many  orders  of  magnitude  of  energy  release  and  are  very  useful  in  an 
empirical  sense.  Scaling  laws  do  not  explicitly  account  for  the  absorption  of  the  air 
or  ground  impedance,  but  instead  simply  state  that  the  distance,  time,  and  impulse 
of  blasts  scale  as  E-173,  or  equivalently  W~1/3.  Here  E  and  W  are  the  energy 
release  and  weight  of  an  explosive,  respectively. 

To  obtain  the  initial  conditions  at  30  m  for  C-4  charge  sizes  other  than  0.57  kg, 
the  blast  durations  can  be  found  from  a  scaling  law.  The  duration  of  the  pulse  from 
a  blast  with  weight  W  should  be  proportional  to  ( W  /W0)1/3  where  W0  is  the  refer¬ 
ence  charge  weight 

Another  relation,73  more  approximate  and  based  on  the  observation  that  blast 
overpressures  decrease  with  distance  at  a  rate  of  about  r-12,  can  be  used  to  predict 


Potential  functions 


t ,  s 


Figure  6.1.  Pulse  spherical  potential  from  Equation  (6.10).  The  solid  line  is  for  the 
exact  I\,  the  dotted  line  is  for  the  approximate  I x. 
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initial  overpressures  for  different  size  blasts  at  a  fixed  distance.  If  the  peak  pressure 

is  known  for  a  charge  of  weight  W  at  a  fixed  range,74  ohe  expression  81og(  W /W0 ) 
gives  the  increase  in  pressure  in  dB  expected  over  that  from  a  charge  of  weight  Wn . 

Using  these  expressions  for  peak  pressure  and  pulse  duration  relative  to  the 
164.1  dB  data.  Table  6.1  gives  values  of  the  approximate  duration  of  the  pulses  r  in 
milliseconds  and  the  weights  W  in  kg,  for  blasts  with  peak  sound  pressure  levels  at 
3  dB  intervals.  The  next  chapters  will  use  these  values  of  peak  SPL  and  r  for  the 
blast  calculation  initial  conditions. 

6.2  Electric  Spark  Initial  Conditions 

Because  of  scaling,  the  small  explosion  relative  durations  given  in  the  last  sec¬ 
tion  are  nearly  the  same  as  those  seen  for  electric  spark  pulses  in  air.  It  turns  out 
that  if  one  divides  the  blast  values  r,  t+ ,  tv  t2,  and  f3  by  1000,  the  wave  shape  is 
similar  to  that  observed  for  an  electric  spark  pulse.  See,  for  example,  the  thesis  of 

Table  6.1.  Scaled  blast  data. 
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Anderson.7''  Here  r  is  approximately  40  (is.  Thus,  the  duration  values  for  sparks 
used  in  Chapter  10  are  r  =  40x  10-6,  tx  =  O.ox  10-6,  £_  =  11.0326x  10~6,  t2  = 
0.50293016x  10“6,  and  t3  =  0.50585915x  10~6  s. 

6.3  Initial  Conditions  near  Surfaces 

Chapters  10  and  11  will  present  propagation  results  including  reflections  from 
hard  surfaces  for  sparks  and  hard  and  porous  surfaces  for  blasts.  This  section 
describes  the  initial  conditions  used  for  those  cases  when  a  reflected  wave  initially 
existed  in  the  computational  domain.  For  simulations  where  an  initial  condition 
contained  a  reflected  wave,  lineai  superposition  was  used  between  the  incident  pulse 
and  its  hard  surface  reflection  image.  This  starting  condition  is  approximate,  but  is 
certainly  satisfactory  for  the  propagation  decay  analysis  in  Chapter  11. 

For  the  high  amplitudes  of  interest  in  this  research,  the  pressure  doubling  seen 
next  to  a  hard  surface  in  linear  acoustics  should  be  replaced  with  a  more  accurate 
relation  dependent  on  the  amplitude  and  past  history  of  the  incident  wave.  Thus 
the  starting  wave  contains  nonphysical  components.  The  nonphysical  wave  com¬ 
ponents,  however,  do  not  severely  affect  the  overall  solution  down  range. 

This  linear  hard  surface  superposition  method  also  was  employed  for  the  initial 
conditions  for  blast  waves  over  a  porous  ground.  Boundary  points  where  the  in¬ 
cident  pulse  initially  intersected  the  ground  were  implemented  using  the  hard 
ground  boundary  condition.  For  all  distances  greater  than  those  initially  touched  by 
the  incident  and  reflected  waves,  the  porous  medium  equations  coupling  the  air  pro¬ 
pagation  equations  acted  as  a  boundary  condition.  See  Figure  6  2.  for  the  initial  po- 
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sitions  in  space  of  the  direct  and  reflected  waves  with  respect  to  the  implementation 
of  the  boundary  conditions  for  this  porous  case. 

The  initial  conditions  described  in  this  chapter  worked  well  for  the  simulations 
of  Chapters  10  and  11.  They  can  be  improved  with  additional  work,  but  as  will  be 
seen  later,  the  conditions  were  satisfactory  for  the  modeling  of  finite  amplitude  pro¬ 
pagation  near  the  hard  and  porous  ground  surfaces. 


Figure  6.2.  The  initial  condition  in  space  for  the  direct  and  reflected  waves  with 
respect  to  the  boundary  conditions  for  the  porous  ground  simulations. 
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7.  NUMERICAL  SOLUTION  TO  THE  AIR  EQUATIONS 

In  the  second  volume  of  his  lectures  on  physics76  Richard  Feynman  states 
“The  only  [sic]  general  methods  of  solution  are  numerical.’’  This  chapter  describes 
the  different  types  of  numerical  solutions  possible  for  the  equations  given  in 
Chapter  4  and  then  specifies  the  particular  solution  chosen  for  this  thesis.  This  type 
of  numerical  solution  should  be  applicable  to  many  types  of  problems  involving 
finite  amplitude  sound  waves  in  fluids,  not  just  the  one  addressed  in  this  work.  The 
properties  of  the  techniques  such  as  accuracy  and  high  frequency  resolution  are  also 
discussed. 

7.1  Background 

A  finite  element  solution  to  (4.39)  is  possible.  However,  an  iterative  solution 
is  required  to  evaluate  the  nonlinear  terms  in  the  system  at  each  time  step.  This 
procedure  is  highly  expensive  from  a  computational  point  of  view.77  New  methods 
of  finite  element  solution  to  nonlinear  problems  are  being  investigated  currently,78 
but  it  seems  that  finite  element  methods  are  best  suited  to  linear  equations. 

Explicit  finite  difference  methods  do  not  have  any  difficulties  with  nonlinear 
equations.  One  can  implement  nonlinear  terms  just  as  easily  as  the  linear  ones  with 
this  technique.  An  explicit  method,  therefore,  will  be  used  in  this  research. 

The  application  of  finite  differences  to  acoustical  equations  is  not  a  recent  idea 
For  example,  an  entire  chapter  is  devoted  to  the  finite  difference  solution  of  sound 
waves  in  Richtmyer  and  Morton’s  Difference  Methods  for  Initial-Value  Problems.'9 
That  work,  however,  did  not  consider  acoustical  nonlinearities.  To  study  blast 
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waves  with  strong  nonlinearities,  Brode  used  a  finite  difference  solution  to  a  set  of 

fluid  dynamics  equations  in  the  1950s.80,81  Brode’ s  methods  were  used  for  overpres¬ 
sures  of  2000  atmospheres  to  0.10  atmosphere  and  are  not  applicable  for  weaker 
waves,  the  subject  of  this  thesis.  Rao  and  Zumwalt  in  1970  modeled  the  pressure 
history  of  a  sonic  boom.82  Their  method  is  not  dissimilar  to  the  one  used  in  this 
thesis.  Rao  and  Zumwalt’ s  finite  difference  method  is  simplistic  and  it  includes  no 
absorption  effects,  but  they  do  solve  a  linearized  system  of  equations.  By  the  early 
1970s  it  was  apparent  that  for  weak  pressure  waves,  the  linearized  compressible  flow 
equations  should  be  used  instead  of  the  full  equations  of  computational  fluid  dynam¬ 
ics.83 

In  1974  Alford  et  al.  stated  that  highly  accurate  fourth-order  finite  difference  al¬ 
gorithms  gave  superior  results  to  those  of  lower-order  for  acoustic  equations.84  This 
“order”  is  meant  in  a  slightly  different  sense  than  the  second-order  nonlinearities 
discussed  earlier.  In  a  finite  difference  approximation,  a  derivative  is  represented  by 
using  a  finite  difference  expression  such  as 


Mx)*  VXx 


where  Ax  is  a  small  increment  in  x.  The  “order”  of  a  finite  difference  approxima¬ 
tion  is  related  to  how  accurately  the  derivative  is  modeled  in  a  numerical  sense. 
Here 


iix  +  Ax )  —  iix) 
Ax 


=  ilfix)  +  Cite  +  C2 Ax2  +  .  .  . 


(7.2) 


where  C  lt  C2,  .  .  .  are  constants.  This  approximation  is  first-order  since  the  approx- 
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imafion  is  dominated  by  the  term  C  jAx .  Another  discretization  of  if(x )  is 


nix  +  Ax )  -  tix  -  Ax ) 


=  n’(x)  +  Z^Ax2  +  D2Ax4  +  .  .  . 


which  is  instead  dominated  by  the  term  Z^Ax2,  so  this  approximation  is  second- 
order.  Similarly  the  approximation 


8iiix  -I-  Ax )  —  8; iix  -  Ax  )  -  iix  -j-  2 Ax )  +  ilix  —  2 Ax ) 

(7.4) 

=  ilf(x)  +E jAx4  +  .  .  . 


is  fourth-order.  Alford  et  al.  also  foimd  that  these  fourth-order  methods  would  al¬ 
low  useful  difference  computations  on  coarse  grids. 

Possibly  unaware  of  the  Alford  et  aL  results,  Walkington  and  Eversman  used 
several  different  second-order  finite  difference  algorithms  to  compute  the  solution  to 
one-dimensional  shocked  acoustic  waves  in  ducts.86,86  They  found  that  their  results 
agreed  with  the  finite  amplitude  acoustic  solution  of  Blackstock  for  an  initially 
sinusoidal  wave.  One  aspect  of  their  work  that  is  particularly  interesting  is  their  use 
of  an  artificial  viscosity  factor.  The  purpose  of  this  artificial  viscosity  is  to  attenuate 
the  high  frequencies  not  well  resolved  by  their  finite  difference  grid.  These  high 
frequencies  are  generated  by  the  finite  amplitude  nonlinear  effects. 

The  most  highly  developed  finite  difference  solution  to  nonlinear  acoustic  equar 
tions  has  been  introduced  by  Maestrello,  Bayliss  and  Turkel.87'91  They  used  fourth- 
order  finite  difference  schemes  to  model  the  two-dimensional  nonlinear  interaction 
of  sound  with  the  shear  layer  of  a  jet  The  particular  technique  they  employed  was  a 
fourth-order  in  space,  second-order  in  time  version  of  the  MacOormack  finite 


47 


difference  method  derived  by  Gottlieb  and  Turkel.92  Because  of  the  completeness  of 
the  work  of  Maestrello,  Bayliss,  and  Turkel,  the  remainder  of  this  dissertation  will 
use  many  of  their  techniques.  S.  I.  Hariharan  also  studied  nonlinear  acoustic  propa¬ 
gation  in  air  proceeding  similarly  to  Maestrello,  Bayliss,  and  Turkel.93  His  methods 
are  highly  mathematical  and  are  not  as  straightforward  as  the  Maestrello,  Bayliss, 
and  Turkel  techniques,  however,  so  the  latter  group’s  methods  will  be  followed. 

7.2  The  Finite  Difference  Method 

For  simplicity  this  research  uses  a  regular  discretization  in  both  the  d  and  z 
directions.  For  any  of  the  variables  the  notation 

ifKt,  d,  z)  =  ipinAt,  iAx,  j tlx)  =  ipp'j  (7.5) 

represents  discretizing  in  time  in  steps  of  At  and  in  space  in  steps  of  Ax . 

The  Maestrello,  Bayliss,  and  Turkel  procedure  begins  by  splitting  Equation 
(4.39), 

wt  -f  Fd  +  Gz  —  SF  +  SG  (7.6) 

into  the  two  equations 

wt  +  Fd  =  SF  (7.7) 

and 

wt  +  Gz  =  SG  .  (7.8) 

One  may  solve  each  split  equation  straightforwardly  by  finite  difference  schemes. 
Call  Ld(At)  and  Lz(At)  finite  difference  solution  operators  which  advances  the 
solution  of  (7.7)  and  (7.8)  one  step,  i.e.,  wit  +  ^f)  =  Ld(At)w(t)  for  (7.7)  and 
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wit  +  At)  =  LziAt)wit)  for  (7.8).  Then  the  Strang  splitting  method  gives  the 
solution  to  ( 7.6)  as 


wit  +  2  At)  =  LdiAt)Lz(At)Lz(At)Ld(At)wit)  (7.9) 


where  both  LdiAt)  and  LziAt)  have  been  applied  twice,  in  a  symmetrical  order.94 
Splitting  is  not  exact,  but  is  merely  another  finite  difference  numerical  approxima¬ 
tion.  The  symmetry  in  (7.9)  is  necessary  to  keep  the  finite  difference  approxima¬ 
tions  high-order. 

The  operator  LA  A t)  is  given  by  two  versions,  the  first  of  which  is 


<  =m/l  +  (-^)(7^-8F/l+1  +  **2)  +&SF? 


<+1  -  j 


wf  +  uf  -SFfLi  +F?_2)  +  AtSF? 


(7.10) 


which  is  a  two-step  method,  fourth-order  accurate  in  Ax  and  second-order  accurate 
in  At .  Here  wn  denotes  a  temporary  value  in  time,  and  Fh  is  the  value  of  F  at  that 
temporary  time.  The  other  v  sion  of  Ld(  At)  is 


w?  =  w?  -(-£L)(7F?  -  +271 2)  +  A tSF? 

bAx 


w?  +  l  =  — 
1  2 


«f  +  uf  +(-^)(7if  -8FR.J  +F?+  2)  +  At  SF? 


(7.11) 


In  using  Ld(  At),  it  is  necessary  to  alternate  the  use  of  (7.10)  and  (7.11)  to  retain  ac¬ 
curacy.96  The  operator  Lz(At)  is  similarly  given  in  two  variants 
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<  -  w?  +  <  -£t)aG?  -8 g;+1  +  o/+2)  +  *  sg; 

<+1  “  I  <  +  -<^r)*7G/  -  8GjL,  +  G/_2)  +  4«SG/ 


(7.12) 


of  =  wf-( -= ~)(7Gf  -  8Gf_1  +  Gf_2 )  +  At  SG/ 


w/+1 


|  “?  +  wf  +  (^-,(7G/  -  80/+1  +  G/+2)  +  *SG/ 


(7.13) 


These  fourth-order  in  Ax  schemes  are  not  applicable  everywhere  on  a  finite 
domain,  however.  The  operator  Ld(  At)  needs  variable  values  at  spatial  points  i— 2, 
£— 1,  i ,  £+1,  and  £+2  to  advance  from  time  n  to  n+1.  If  £  runs  from  1  to  7,  the 
above  routines  are  not  applicable  at  points  1,  2,  7—1,  and  7.  The  points  1  and  7  are 
boundary  points,  and  numerical  boundary  conditions  given  in  Chapter  8  will  ad¬ 
vance  their  values  from  time  n  to  time  n+1. 

For  the  variable  values  at  points  i  =  2  and  i  =1—1  the  second-order  in  At  and 
Ax  versions  of  Ld(At)  were  used  in  this  research  Just  as  with  the  fourth-order  in 
Ax  schemes,  the  second-order  Lrf( At)  has  two  variants 


<  -  «f  +  ~  *T+1)  +  &SF? 

w?+l  =  ~  +  +  At  SFf 

4  OX 


(7.14) 


wf  =  wp~  -=—(Fp  -  TJLi )  +&SF? 
OX 


v?+1  =  ~  wf  +  w?  +  J*4Ff-Ff+l)  +  AtSTf 

4  OX 


(7.15) 
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which  one  should  alternate.  Similarly  at  j  =  2  and  j  —  J—l  for  1  <  j  <  J,  the 
second-order  L2(  A )  schemes  are 


f  -  wf  + -%AG?  -  G?+l)  +  X  SG? 


■i  =  l  w*  +  wf  -  -=1( Gf  -  Gf_x )  +  A  SGf 
2  J  J  Ac  J  J  J 


(7.16) 


wf  =  wf  -  ~{Gf  -  Gf_i )  +  A  SGf 

wf+1  =  \wf  +  wf  +  jtiGf  -  Gf+l )  +  ASG?  . 


(7.17) 


In  the  future  it  might  be  possible  to  find  some  higher-order  schemes  which  apply  at 
i  —  2  and  i  =1—1  and  j  —2  and  j  =  J—l,  but  the  second-order  schemes  work 
well. 


7.3  Phase  Dispersion  and  Artificial  Viscosity 

Note  the  reason  that  the  fourth-order  schemes  are  better  for  wave  propagation 
than  the  second-order  schemes  is  because  the  former  have  minimal  phase  disper¬ 
sion.  The  one  way  model  equation  ~  +  c0  ~  =  0  will  serve  as  an  example.  If 

at  ox 

one  substitutes  the  harmonic  solution p(t,x)  =  exp[r(A  — fee)],  the  model  equation 
reduces  to  the  dispersion  relation  w  =  c0k. 

However,  if  one  substitutes  the  same  time  harmonic  solution  into  a  finite 
difference  approximation  to  the  model  equation,  a  dispersion  relation  w  =  c0k  +  ET 
results  where  ET  are  error  terms.  For  example  the  second-order  in  Ac  and  A 
MacCormack  method  applied  to  the  model  equation  gives 
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pr"  =  p?  -  cf^n\ i 

1  co& 

+  2"  A^-  (p/+ 1  “  +  Pi- 1  )  . 


(7.18) 


which  turns  out  to  be  a  well  known  finite  difference  method,  the  Lax-Wendrofif 
scheme.  If  the  substitution  pilAx,  nA£)  =  expti(^nA£  —  klAc)]  is  made,  the  ex¬ 
act  dispersion  relation  for  second-order  MacCormack  is 


Y2 

—i(elujU  —  1)  =  c0  -^sin(&  Ax)  +2 i  cQ-^~  sin2(&Ax/2) 
Ax  0  Ax 


(7.19) 


Similarly  the  fourth-order  MacCormack  method  applied  to  the  model  equation  gives 


pT1  =Pl  +  2  -  Sp/Vl  +  *Pl- 1  -Pi-2) 


_1_  ^ 
72  C°  Ax 


(7.20) 


(7p/*+2  -  64p/Vi  +  11 4pf  -  64p/Li  +  7p/L2) 


which  has  an  exact  dispersion  relation 


—i(el“u  -  1)  = 


C°  ^  8sin(fc  Ax)  —  sin(2&  Ar )  1 
6  Ax  J 


.  f  |!f  i 

±-  Co-=-  64sin2(*  Ax/2)  -  7an2(kAx) 

18  0  Ax  «■  J 


(7.21) 


The  derivation  is  straightforward.  Notice  that  both  Equations  (7.19)  and  (7.21)  ap¬ 


proach  ijj  —  c0k  as  At  and  Ax  approach  zero. 


These  dispersion  relation  expressions  were  expanded  in  Taylor  series  for 
sin(&  Ax)  and  exp(ic*;Af]  with  the  symbolic  manipulation  program,  Mathematica96 
Using  these  Taylor  series  expansions,  the  expressions  for  the  error  terms  ET  in  the 
dispersion  relations  w  =  c0k  +  ET  were  found  for  both  the  second-order  and 
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fourth-order-  finite  difference  schemes.  Figure  7.1  plots  the  relative  magnitude  of 
ET  for  the  two  cases  using  constants  appropriate  for  electric  spark  pulses,  c0  =  343 
m/s,  Ac  =  250x  10~6  s  ,  and  A/  =  200x  10~9  s.  Note  that  the  values  of  | ET  |  for 
the  fourth-order  scheme  are  less  than  those  for  the  second-order  scheme  for  all 
wavenumbers,  and  that  the  difference  between  the  two  errors  increases  with  increas¬ 
ing  wavenumber.  The  idea  for  plotting  the  errors  of  dispersion  relations  was  taken 
from  the  dispersion  relation  plots  in  the  Ph.  D.  thesis  of  Trefethen.97,98 

The  improvement  from  the  second  to  the  fourth-order  methods  is  possible  only 
because  the  ratio  c0  At /Ac  is  nearly  0.25.  If  one  uses  much  higher  or  lower  values 
of  c0AT/Ar,  the  fourth-order  method  will  lose  accuracy  since  it  is  fourth-order  in 
Ac  but  only  second-order  in  At.  By  keeping  c0  At  small  in  a  relative  sense  to  Ac, 
the  errors  a°sc elated  with  At2  will  not  dominate  the  errors  associated  with  Ac4.  See 
the  Turkel  reference  for  the  details  on  this  criterion. 


Even  with  the  fourth-order  in  Ac  MacCormsck  method,  the  nonlinearities  in 
the  finite  amplitude  equations  produce  high  frequencies  that  the  finite  difference 
methods  cannot  well  resolve.  Such  inaccuracies  will  be  evident  near  discontinuities 
such  as  a  shock  front,  where  the  variables  will  change  quickly  over  a  few  spatial 
steps.  High  frequency  dispersion  can  be  minimized  by  the  use  of  an  artificial  viscos¬ 
ity.  As  mentioned  earlier  in  this  chapter,  Walkington  and  Eversman  have  used  a 
fourth-order  artificial  viscosity.  Their  method  takes  the  form 


{/•/*  +  ! 

^  l 


|new 


i/i»  +  1 


old  ~  U  ^*+2  -  +  6^f  -  4^f-l  +  ^"-2 


( 7.22) 


where  i>  is  small.  One  way  of  thinking  about  this  process  is  that  a  low  pass  filter  is 
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1000.  2000.  3000.  4000.  5000.  6000. 


Figure  7.1.  The  relative  magnitudes  of  the  error  terms  ET  versus  wavenumber  for 
two  finite  difference  schemes.  The  second-order  MacConnack  method  is 
represented  by  the  solid  line,  and  the  fourth-order  in  space  and  second-order  in 
time  MacCormack  method  is  given  by  the  dashed  line. 
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run  on  the  variable  d  after  each  time  step.  More  correctly,  note  that  the  above 

effect  is  equivalent  to  adding  an  artificial  viscosity  to  the  model  equation 

dll'  dll’  A  ,  . 

- 1-  cn  — —  =  0  to  obtam 

dt  d X 


One  usually  incorporates  classical  shear  viscosity  into  the  model  equation  as 


dtl'  dil' 

~dT  +  c°  aT  “  • 


(7.25) 


Ik  is  well  known  that  absorption  in  (7.25)  is  proportional  to  uA  Similarly  one  can 
show  with  the  assumptions  ^  close  to  c0k  and  /iartu.’3/c04<<  1,  that 
exp[i  ivJt  —  kx )  ]  substituted  into  ( 7.22)  yields  k  =  k0  —i  where 
a  an  =  /iartU.'Vco6.  Hence,  adding  a  fourth-order  artificial  viscosity  adds  absorption 
which  is  proportional  to  u/4.  By  picking  carefully,  one  should  be  able  to  dissi¬ 
pate  the  high  frequencies  not  well  resolved  by  the  spatial  grid,  while  leaving  the 
lower  frequencies  nearly  unaffected. 


To  implement  the  viscosity  technique,  Equation  (7.22)  was  used  on  those 
points  in  the  computational  domain  that  were  calculated  using  the  fourth-order  in 
Ja  MacCormack  scheme.  The  method  was  employed  every  fime  step  in  the  follow¬ 
ing  four  places:  in  the  first  and  second  components  of  Equation  (7.7)  and  in  the  first 
and  third  components  of  Equation  ( 7.8).  This  procedure  applies  the  artificial  viscos- 
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ity  to  both  of  the  split  components  of  the  continuity  equation  and  to  each  of  the 
particle  velocity  equations.  Chapter  10  will  provide  some  specific  values  of  v  that 
were  used  during  calculations. 

7.4  Dissipation,  High  Frequency  Resolution,  and  Implementation 

Another  advantage  of  using  the  MacCormack  finite  difference  schemes  is  that 
these  methods  inherently  contain  their  own  dissipative  terms,  which  under  usual  cir¬ 
cumstances  are  completely  ignored.  This  dissipation  is  extremely  useful  in  killing 
parasitic,  nonphysical  waves  that  can  be  launched  by  the  interface  of  a  difference 
scheme  and  any  boundary  condition,  which  can  be  absorbing,  hard  or  porous.  Such 
waves  often  take  the  form  of  (— l)n,  (—1)*,  and  (—  IV  and  can  quickly  swamp  a 
valid  numerical  solution  by  exponential  growth-  Trefethen  investigated  such  waves 
in  detail.  Sparrow’s  previous  attempt  to  solve  acoustical  equations  was  beset  by 
problems  of  this  nature.  That  investigation  used  the  Leapfrog  finite  difference 
scheme,  a  simple  and  standard  numerical  method."  But  Leapfrog  is  inherently  non- 
dissipative,  which  renders  it  useless  for  realistic  computations  involving  absorbing 
and  surface  boundary  conditions. 

In  using  the  high-order  in  space  finite  difference  schemes  of  the  preceding  para¬ 
graphs,  it  may  not  be  clear  how  well  the  high  frequencies  will  be  resolved.  As  a 
rule  of  thumb,  one  should  always  use  at  least  10  to  12  spatial  grid  points  for  one 
wavelength  of  the  highest  frequency.  Once  the  spatial  grid  size  Ax  has  been  set  us¬ 
ing  this  rule,  the  time  step  AT  is  then  known  from  the  necessity  of  keeping  c0  Ai/Ax 
at  about  0.25,  as  was  discussed  previously. 
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Now  examine  the  requirements  for  resolving  the  high  frequencies  involved  in 
the  initial  near  discontinuity  of  a  blast  pulse.  Suppose  a  blast  wave  contains  a  near 
discontinuity  which  has  a  rise  time  of  rt  =  lx  10-4  s,  with  frequencies  up  to  10 
times  the  1  Art)  frequency.  Frequency  resolution  up  to  100  kHz  is  required  to 
model  this  wave.  For  a  sound  speed  of  343  m/s,  the  wavelength  of  100  kHz  is 
3.43x  10~3  m,  and  Ac  should  be  no  larger  than  3.43x  10~4  m. 

In  the  results  of  Chapters  10  and  li,  the  high  frequencies  contained  in  the 
spark  and  blast  pulses  will  not  be  well  resolved  due  to  computer  memory  limita¬ 
tions,  so  the  solutions  will  not  be  accurate  in  the  vicinity  of  the  near  discontinuities. 
With  the  present  computer  program  one  cannot  practically  perform  numerous  high 
frequency  rise  time  calculations  in  two  dimensions  over  large  spatial  grids. 

The  above  algorithms  have  been  incorporated  into  a  highly  vectorized  computer 
program  in  Fortran.  The  rectangular  symmetry  of  the  spatial  grid  and  the  two  step 
fourth-order  difference  schemes  are  ideal  for  parallel  implementation  on  modem 
supercomputers.  Computationally  this  program  requires  up  to  4  to  5  megawords 
{  million  floating  point  numbers)  of  continuous  memory  usage  and  up  to  20  minutes 
of  CPU  time  on  a  CRAY-2  for  medium  sized  runs.  If  one  halves  the  spatial  grid 
size  to  better  resolve  high  frequencies,  one  needs  4  times  the  memory,  since  this  is 
spatially  a  2-D  program,  and  2x4  =  8  times  the  CPU  time.  If  one  is  interested  in 
propagation  primarily  in  the  d  direction,  the  program  can  window  a  pulse  moving  in 
that  direction.  This  feature  has  already  cut  the  memory  usage  of  the  simulation 
drastically.  However,  adaptive  grids  or  different  difference  schemes  may  allow  one 
to  reduce  this  memory  burden  in  the  future. 
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Another  important  point  is  that  in  a  typical  run,  the  program  will  generate 
literally  millions  of  numbers.  Public  domain  software  written  at  the  National  Center 
for  Supercomputing  Applications  (NCSA)  was  used  to  translate  these  numbers  into 
color  images  for  interpretation  on  videotape  and  color  graphics  terminals.100  In  addi¬ 
tion,  the  program  allows  the  user  to  place  “numerical  receivers”  anywhere 
throughout  the  spatial  domain.  These  receivers  can  record  the  pressure  waveforms 
at  the  specified  spatial  points  for  the  duration  of  a  run.  The  present  code  is  over 
5000  lines  long,  and  Appendix  B  describes  it  in  detail. 
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8.  NUMERICAL  ABSORBING  BOUNDARY  CONDITIONS 

This  chapter  describes  the  numerical  absorbing  boundary  conditions  that  are 
used  in  this  dissertation.  First  some  background  on  these  conditions  is  given.  The 
radiation  conditions  used  in  this  research  to  correctly  model  outgoing  waves  are 
then  detailed.  Where  radiation  conditions  are  inappropriate,  characteristic  condi¬ 
tions  are  used,  and  these  are  described  also.  One  might  use  the  characteristic  condi¬ 
tions  for  radiation  conditions,  but  it  is  shown  that  the  ones  first  given  in  this  chapter 
are  superior. 

Figure  8.1  shows  the  computational  domain  used  in  this  thesis.  It  is  bounded 
by  the  four  sides  z  =  zmm ,  z  =  zmax ,  d  ~  dmin,  and  d  =  cfmax  in  the  (d,z)  plane. 
The  source  of  sound  is  assumed  to  be  at  the  origin  id  =  0,2  =0).  For  free  field 
propagation,  radiation  conditions  are  used  for  all  of  the  boundaries  except  ci  =  dmw 
where  the  characteristic  condition  is  used.  For  the  blast  modeling  of  Chapter  11, 
the  boundary  condition  along  z  =  2min  is  changed  from  radiation  to  either  that  for  a 
hard  or  a  porous  ground  surface. 

8.1  Development  of  Absorbing  Boundary  Conditions 

First  of  all,  what  are  absorbing,  or  radiation,  boundary  conditions?  The  solu¬ 
tion  to  acoustical  equations  using  finite  differences  is  restricted  to  a  finite  domain 
since  computers  have  limited  memory.  Often  a  larger  domain  is  modeled,  however, 
and  boundary  conditions  on  the  finite  computational  domain  must  be  used  which  al¬ 
low  waves  to  pass  out  of  the  computational  domain  without  reflections. 
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The  simulation  of  nearly  reflectionless  boundaries  has  been  of  recent  interest  in 
many  fields.  The  two  approaches  used  are  the  placement  of  an  absorbing  region 
surrounding  the  computational  domain  to  slowly  “eat  up”  the  waves  which  impinge 
on  it  or  the  use  of  an  absorbing  or  radiation  condition  which  is  local  in  space  at  the 
edge  of  the  domain  and  which  exactly  satisfies  an  outgoing  condition.  The  approach 
taken  in  this  thesis  focuses  on  this  latter  method  since  the  former  requires  more 
computer  memory. 

The  problem  of  false  reflections  from  absorbing  boundaries  has  been  overcome 
in  recent  years  for  scalar  plane  wave  propagation.  Researchers  have  been  able  to 
develop  nearly  reflectionless  conditions  for  model  equations  such  as 

P'tt  +co2(PrE  +  P'yy^  (8.1) 

for  plane  waves  with  arbitrary  angles  of  incidence  on  the  boundaiy.  Here  again  p'  is 
the  acoustic  density  deviation.  Higdon101,102  has  proposed  difference  approxima¬ 
tions  to  the  boundaiy  condition 


ft 

7-1 


,  d  d  .  ,  n 

(cos a.  —  -  c0  —)p  =  0 
at  etc 


(8.2) 


which  is  reflectionless  for  the  2 p  plane  waves  at  angles  ±ocj.  Smilar,  less 
developed  conditions  have  been  given  by  Keys.103  Higdon  has  shown  that  the  often 
cited  boundaiy  conditions  of  Engquist  and  Mqjda104,106  are  actually  those  in  (8.2) 
with  =  0  for  all  j ,  i.e.,  they  exactly  satisfy  plane  waves  only  at  normal  incidence. 


These  results,  however,  are  not  sufficient  to  find  an  exact  condition  for  a  spher¬ 


ical  wave  propagating  in  the  system  of  equations 
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+  Po  ( ud  +  vz  +  -  0 

(8.3) 

Po^t  +c02p'd  =  0 

(8.4) 

Po  vt  +  C0Vz  =  0  . 

(8.5) 

These  are  the  resulting  two-dimensional  linear  lossless  equations  presented  in 
Chapter  4  with  the  equation  of  state  substituted  into  (8.4)  and  (8.5).  Gottlieb  et 
al.106  have  demonstrated  that  the  scalar  results  can  be  used  in  the  case  of  a  hyper¬ 
bolic  system  if  the  boundary  conditions  are  applied  to  the  characteristic  variables  of 
the  system,  instead  of  to  the  natural  variables  such  as  p,  u,  and  v  in  (8.3),  (8.4), 
and  (8.5).  Characteristic  variables  are  those  variables  which  decouple,  or  nearly 
decouple,  the  equations  in  a  system  of  equations  from  each  other.  Coughran107  also 
has  stated  that  one  should  base  boundary  conditions  on  characteristic  variables 
whenever  possible.  But  the  Gottlieb  et  al.  result  which  extends  the  conditions  for 
scalar  equations  to  systems  of  equations  yields  reflectionless  boundary  conditions 
only  for  the  case  of  plane  waves,  not  spherical  waves.  Additional  insight  into  the 
complex  problem  present  in  finding  boundary  conditions  for  hyperbolic  systems  can 
be  gained  from  Higdon.108 

Bayliss  and  Turkel  have  examined  absorbing  boundary  conditions  for  spherical 
waves  in  a  hyperbolic  system  of  acoustics  equations.109,110  They  propose  using  a  se¬ 
quence  of  boundary  operators  Bn  to  successively  annihilate  the  first  n  terms  of  the 
expansion 

p(t,r,6,<t>)  = 

j- 1 


(8.6) 
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The  Bx  operator  is  useful  for  the  purpose  of  this  thesis  and  will  be  used  later.  In 
addition,  Hagstrom  and  Hariharan  have  developed  spherical  wave  absorbing  boun¬ 
dary  conditions  for  the  Euler  equations,111  but  these  are  not  useful  for  the  particular 
problem  addressed  here. 

Absorbing  boundary  conditions  must  be  stable  numerically  as  well  as  accurate, 
but  finding  stable  conditions  is  not  straightforward.  The  recent  work  of  Gustafs- 
son112  for  example,  shows  that  high-order  boundary  conditions  lead  to  ill-posed 
problems  (unstable  numerical  schemes).  Therefore,  only  conditions  of  order  one 
have  been  used  in  this  research. 


8.2  Radiation  Boundary  Conditions 


To  exactly  absorb  a  pulse  from  a  monopole  spherical  source  in  Equations  ( 8.3) , 
(8.4),  and  (8.5),  notice  that  the  acoustic  density  p'  exactly  satisfies  the  separable 
wave  equation 


irp’)  =  0 


(8.7) 


while  u,  v,  and  vr  =  usin#  -f-  vcosQ  do  not  Here,  vr  is  the  spherical  component  of 
particle  velocity  in  the  spherical  radial  direction  r.  Because  the  particle  velocities 
will  not  satisfy  (8.7),  the  conditions  on  the  density  and  particle  velocities  will  be 
different  Taking  the  outgoing  part  of  (8.7)  gives  the  condition 


d_ 
d t 


+  l> 

r 


P'  =  0 


(8.8) 


which  is  the  Bayliss  Turkel  Bx  condition  on  the  acoustic  density.  Now  look  at  the 
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one-dimensional  equations  similar  to  (8.3),  (8.4),  and  (8.5)  in  the  spherical  coordi¬ 
nate  system  ir,<t>,9)  with  only  r  variations. 


2vr  dvr 

r  ^  dr 


+  Po  ——  +  -5T-  =  0 


(8.9) 


dt  0  dr 


(8.10) 


Equation  (8.10)  gives  the  exact  condition  relating  p'  to  the  time  derivative  of  vr. 

It  is  not  safe  to  base  a  boundary  condition  on  (8.10)  explicitly,  because  it  in¬ 
volves  vr  which  depends  on  both  u  and  v.  For  boundaries  where  z  or  d  are  con¬ 
stant,  as  in  Figure  8.1,  the  boundary  conditions  should  not  not  involve  the  variables 
u  and  v,  respectively.  This  is  one  of  the  results  of  an  important  paper  by  Mqjda  and 
Osher,113  and  it  is  important  for  this  thesis  since  boundaries  of  the  type  z  or  d  a 
constant  will  be  used. 

For  an  upper  boundary  z  =  zmax  with  the  interior  domain  being  zmax  <  z,  one 
can  insert  the  geometrical  relation  v  =  vrcosd  in  (8.10)  to  produce  the  equation 


uu  ,  n  ~o  uu 

+  cose - £- 

dt  P0  dr 


(8.11) 


C\  n  f 

The  term  is  easily  found  from 
dr 


dp'  Q  dp'  .  a  dp' 

~~  =  cos $-£-  + 
dr  dz  da 


(8.12) 


Thus,  valid  absorbing  boundary  conditions  for  z  =zmax  are  (8.8)  and  (8.11) 
with  the  insertion  of  (8.12) .  They  are  stable  and  exact  for  a  mono  pole  source  at  the 
origin.  This  condition  is  also  correct  for  a  boundary  z  =  z^ n,  except  that  its  imple- 
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mentation  is  different.  The  implementation  used  for  zmax  absorbing  boundary  con¬ 
dition  was 


,,n  +  l  ,,n 
viJ  ~vi,J 


COS  6- 


»  2 
0 

Po 


+ 


n'n  rt'/I  n'K 

Pi.J-Pi-l,J  .  Q  ,  Pi,J~Pi,J- 1  Q 

— sm 6  -| - - - cosy 


Id 


(8.13) 


=  0 


for  2<j<I  in  the  d  direction  in  the  domain  1  <i</  in  d  and  1  <j<J  in  z.  For 
2min  the  implementation  was 


,,n  +  l  _  ,,n 
vi, 1  vi, 1 


+ 


cosd- 


Po  L 


n’n  _  r,’n  n  __  „'n 

P  i,l  P  i-1,1  .a  P  i,2  P  i,l  „ 
—sin#  + - 7 - —cosd 


/\z 


(8.14) 


=  0 


In  the  case  of  an  outgoing  condition  at  d  =  dmax,  where  u  is  needed  instead  of 
v ,  the  relation 

~  +  sin#— ^  «  0  (8.15) 

dt  p0  dr 

was  similarly  obtained  from  (8.10).  Equations  (8.15)  and  (8.8)  were  used  for  the 
absorbing  boundary  here.  The  particle  velocity  condition  for  d  =  dmax  was  imple¬ 
mented  as 


“/71  -  vlj 


& 


+ 


sin#  — 
Po  L 


^  "  ^'L^o  +  PliZjljzl  cq^ 


dd 


(8.16) 


=  0  . 
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For  each  of  the  boundaries  z  =  zmax,  z  =  2min*  and  d  =  dmax,  Equation  (8.8)  was 
discretized  similarly. 

All  of  the  boundaries  derived  in  this  chapter  are  for  linear  lossless  conditions. 
Fortunately,  the  absorption  and  nonlinearities  present  in  the  air  propagation  equar 
tions  are  weak.  This  fact  allows  the  use  of  the  linear  conditions  specified  here  to 
work  with  the  nonlinear  lossy  equations  solved  in  this  thesis. 

The  conditions  here  do  not  account  for  the  variable  Sfr ,  the  acoustic  entropy 
deviation,  however.  But  since  this  variable  is  associated  only  with  the  losses  in  the 
equations  of  Chapter  4  and  7,  its  specification  is  not  critical  near  absorbing  boun¬ 
daries,  but  is  so  near  physical  surfaces.  Near  the  boundary  d  =  dmin  and  near  the 
radiation  boundaries,  Sfr  was  updated  by  simple  extrapolations.  Chapter  10  will 
specify  boundary  conditions  on  Sfr,  as  well  as  p\  u,  and  v  for  hard  surfaces. 

8.3  Characteristic  Boundary  Conditions 

The  boundary  at  d  —  dmiD  is  not  a  radiation  boundary,  since  in  this  thesis  the 
sound  pulses  will  be  initially  placed  at  d  >  dmiD  and  propagate  in  the  -k/  direction. 
The  boundary  condition  implemented  along  d  —  dmin  is  called  a  characteristic  boun¬ 
dary  since  it  is  based  on  the  concept  of  characteristic  variables,  defined  in  Section 
8.1. 

Recall  the  one-dimensional  spherically  symmetric  acoustic  equations  in  density 
deviation  and  radial  particle  velocity,  Equations  (8.9)  and  (8.10).  It  is  possible  to 
rewrite  iht-ac  cvjiuuions  in  the  form 
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w'f  + 


c0  0 
0  — c 


co  1  ~ 1 

<  +  —  i  _i  w'  * 

o  r  -*•  x 


where 


(8.17) 


W  = 


1  ,  .  Po 

— —p  +  - 

I V2  c„V§  r 


v5p 


( 8.18) 


and  where  the  coefficient  matrix  of  v/r  is  diagonal. 

This  transformation  is  made  possible  by  writing  (8.9)  and  (8.10)  as  the  system 


Af  +  BAf.  +  CA  =  0 


(8.19) 


where 


0  Po 
C  ^ 

^  o 


(8.20) 


(8.21) 


(8.22) 


.0  0  • 


The  eigenvalues  of  B  are  easily  found  as  c0  and  — c0 ,  and  two  right  eigenvectors  as¬ 
sociated  with  these  eigenvalues  are 


r\  - 


(8.23) 


and 
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r2  = 


(8.24) 


respectively.  One  can  also  find  two  left  eigenvalues  associated  with  c0  and  —c0 . 
They  are 


,  _  1  Po 

h~  ^ 


(8.25) 


1  -Po 
v5  c0v5  ’ 


(8.26) 


respectively.  Finding  such  eigenvectors  means  B  may  be  diagonalized  through 


A  =SBS“ 


(8.27) 


where 


(8.28) 


r2] 


(8.29) 


Multiplying  (8.19)  by  S  from  the  left  and  expanding  A  as  A  —  S_1SA  =  S_1W  gives 
( 8.17).  Here  the  variables  W  are  related  to  A  by 


W  =  SA 


(8.30) 


A  =  S“  V  . 


(8.31) 
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The  variable  wV  is  called  the  outgoing  characteristic  variable  of  these  equations,  and 
w 2'  is  the  incoming  characteristic  variable. 

Similarly,  one  can  diagonalize  the  cylindrical  2-D  Equations  (8.3),  (8.4),  and 
(8.5),  with  respect  to  the  d  direction.  The  resulting  equations  are 


w,  +  Dwj  +  EjWz  +  E2w  =  0 


for  a  diagonal  D  matrix  where 


and 
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(8.32) 


(8.33) 


(8.34) 


(8.35) 


(8.36) 


Here  the  characteristic  variable  with  respect  to  the  +d  direction  is  w  1,  the  charac- 
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teristic  variable  with  respect  to  the  —d  direction  is  w 2,  and  w 3  =  v  is  not  a  charac¬ 
teristic  variable.  Examine  (8.32)  for  the  case  of  plane  wave  propagation  in  the  +  or 
—  direction.  Here  v  =  w 3  =  0.  For  propagation  in  the  +d  direction,  the  density 

Po 

and  particle  velocity  are  in  phase  as  p'  =  — u ,  which  means  w 2  =  0.  Similarly,  for 

C0 

the  —d  direction,  p'  —  ——u,  so  wl  =  0. 

Co 

Because  of  these  properties,  the  second  case  wl  —  0  was  used  as  the  boundary 
condition  nearest  the  source  at  d  =  dmin  >  0  for  this  research.  Specifying  ml  =  0 
implies  that  no  waves  will  enter  the  domain  d  >  dmin  in  the  direction  +d.  To  im¬ 
plement  holding  wl  —  0,  the  variable  w2  was  extrapolated  over  one  in  the  d 
direction  along  the  boundary  d  =  dmm.  In  finite  difference  notation 

m2{l+1  =  w2lj  (8.37) 

for  all  j  implements  the  method,  using  the  first  two  components  of  (8.36)  and  the 
first  two  components  of 

^-(u;l  +  w2) 

— %=~(ml  —  a;  2) 
w3 

with  wl  =  0.  The  variable  v  is  not  changed  on  a  d^  boundary  update. 

Now  return  to  the  1-D  spherical  Equations  (8.17)  and  (8.18).  Unlike  the  plane 
wave  case,  p'  and  v  are  out  of  phase  except  in  the  very  far  field.  Restated,  an  out¬ 
going  spherical  wave  contains  both  an  incoming  and  an  outgoing  characteristic  com- 


(8.38) 


70 


ponent  for  finite  distances  from  the  source.  This  fact  means  that  if  one  sets  u>2'  =  0 
in  (8.17)  and  (8.18),  that  the  first  equation  of  (8.17)  used  as  an  absorbing  boundary 
condition  will  produce  false  reflections  for  outgoing  waves.  To  demonstrate  this 
point,  the  next  section  will  contain  a  comparison  between  the  condition  w 2'  =  0  as  a 
radiation  boundary  condition  and  the  conditions  derived  in  Section  8.2. 

8.4  Absorbing  Boundary  Condition  Comparison 

This  section  compares  three  candidates  for  absorbing  boundary  conditions  to 
see  which  has  the  minimum  reflection.  The  upper  radiation  boundary  condition, 
2  =  2 max*  was  the  boundary  tested. 

The  equations  of  Chapter  4  were  solved  using  the  methods  of  Chapter  7  in  this 
comparison,  under  the  stipulations  k,  p,  and  pg  =  0,  with  a  peak  SPL  of  140  dB. 
This  is  equivalent  to  solving  a  linear  lossless  system  of  acoustics  equations.  The 
upper  boundary  was  placed  at  zmax  =  0.03  m,  and  the  computational  domain  was 

2  ~  2 max- 

A  spherical  spark  pulse,  of  the  type  described  in  Chapter  6,  was  used  as  an  ini¬ 
tial  condition  on  p\  u,  and  v.  Initially  the  pulse  extended  to  0.03  m  from  the  ori¬ 
gin,  with  c0  =  343  m/sec  and  p0  =  1.21  kg/m3,  values  appropriate  for  air. 

Four  computations  were  performed  with  different  conditions  used  along 
2  ~  2 max-  hi  311  four  calculations  values  of  u  were  found  from  the  interior  equa¬ 
tions,  since  this  is  not  a  boundary  variable  along  z  —  zmax.  Computations  one 
through  four  varied  in  the  way  p'  and  u,  the  boundary  variables,  were  calculated. 
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In  case  one  the  spherical  characteristic  variables  of  Section  8.3, 

dwl'  dwY  ,  Cc  (  0/)  „  /oon> 

—  +  c0  — - (- — \wY-w2\  =  Q  (8.39) 

at  cjr  r  t  ' 

and  w2'  =  0,  were  used.  Case  two  used  the  Bayliss  Turkel  Bx  boundary  condition 
on  the  density  deviation, 

p't  +c0(p'r  +  p’/r)  =  0  ,  (8.40) 

and  the  interior  equation  on  v, 

P0vt+c02p'z=  0.  •  (8.41) 

Case  three  also  used  the  Bl  condition  on  p\  but  on  v  used 

p0vt  +  cos 0c2p'r  =  0  ,  (8.42) 

the  pair  of  boundary  conditions  from  Section  8.2.  The  fourth  case,  to  act  as  the 
ideal  condition,  was  implemented  by  extending  the  domain  beyond  2max  in  z,  so 
that  there  would  be  no  upper  boundary  to  reflect  off.  The  computation  time  for 
case  four  was  significantly  greater  than  that  of  the  other  three  cases  because  of  the 
former’s  larger  domain.  The  conditions  in  cases  one  through  three  were  discretized 
using  simple  forward  and  backward  differences. 

For  the  four  computations,  Sx  =  2d  =  2z  =  250xl0-6  m  and  = 
200*  10-9  s  were  used  making  c0  2t/ Ar  =  0.2744  which  is  appropriate  for  the 
fourth-order  method.  The  pulse  was  allowed  to  propagate  for  600  time  steps  or  120 
microseconds  from  its  initial  condition.  A  numerical  receiver  was  placed  at  (d  = 
0.04  in,  z  =  0.025  m)  to  monitor  the  wave  passage.  Thus  the  receiver  intercepts 
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the  reflected  part  of  the  wave  incident  on  the  absorbing  boundary  at  about  50  de¬ 
grees  from  normal.  The  computation  was  performed  on  a  Cray-2  supercomputer, 
and  each  of  the  four  runs  took  about  180  s  on  a  single  CPU. 

Figures  8.2,  8.3,  and  8.4  are  plots  of  acoustic  pressure,  p  =  c02p',  at  the  re¬ 
ceiver,  comparing  the  ideal  condition  with  the  boundary  conditions  one  through 
three,  respectively.  Case  four  is  plotted  in  each  of  the  three  figures  with  a  dashed 
line,  and  each  of  the  others  is  represented  by  a  solid  line.  Pressure  is  given  in  Pas¬ 
cals  and  time  in  microseconds.  Figures  8.5,  8.6  and  8.7  show  the  difference 
between  conditions  one  through  three  and  the  ideal  condition  four.  For  example. 
Figure  8.2  shows  condition  one  with  a  solid  line,  the  ideal  condition  four  with  a 
dashed  line,  and  Figure  8.5  gives  their  difference. 

Etch  of  the  computations  was  identical  until  about  65  ps  when  the  wave 
reflected  from  the  artificial  boundaries  arrived  at  the  receiver.  Figures  8.2  and  8.5 
show  the  reflection  over  much  of  the  trough  of  the  pulse  for  condition  one,  which  is 
due  to  the  continual  misspecification  of  the  incoming  characteristic  variable,  =  0. 
Hence,  this  comparison  verifies  that  in  the  near  field  there  is  an  incoming  charac¬ 
teristic  component  for  an  outgoing  spherical  wave. 


Figure  8.2.  Acoustic  pressure  in  Pa  versus  time  at  the  numerical  receiver, 
one  is  the  solid  line,  and  case  four,  the  ideal,  is  dashed. 
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Figure  8.3.  Acoustic  pressure  in  Pa  versus  time  at  the  numerical  receiver.  Case  two 
is  the  solid  line,  and  case  four,  the  ideal,  is  dashed. 
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Figure  8.7.  Acoustic  pressure  difference  in  Pa  between  cases  three  and 
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Case  two  in  Figures  8.3  and  8.6  shows  a  substantial  deviation  from  the  ideal 
which  is  greatest  at  the  beginning  of  the  reflection,  while  case  three,  the  conditions 
of  Section  8.2,  is  nearly  identical  to  the  ideal  in  Figures  8.4  and  8.7.  The  dips  in  the 
waveforms  in  front  and  behind  the  peak  of  the  pulse  for  all  four  cases  are  not  due 
to  reflections.  They  were  due  solely  to  numerical  dispersion  in  the  interior  finite 
difference  method.  It  is  near  large  waveform  gradients  that  the  high  frequencies  not 
being  well  resolved  by  the  grid  first  manifest  themselves. 

Not  shown  in  Figures  8.2  through  8.7  is  the  fact  that  one  could  see  almost  no 
reflected  waves  off  the  numerical  boundary  for  case  three  in  color  raster  plots  of  the 
numerical  domain.  On  the  other  hand,  the  case  one  boundary  always  showed  a 
strong  spherical  wave  reflection,  and  the  case  two  boundary  showed  a  stronger 
reflection  the  more  grazing  the  angle  of  incidence  was,  although  for  near  normal  in¬ 
cidence  little  reflection  was  present 
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9.  NUMERICAL  SOLUTION  TO  THE  POROUS  EQUATIONS 

This  chapter  will  present  a  numerical  solution  for  the  time-dependent 
phenomenological  equations  for  the  porous  medium  described  in  Chapter  5.  This 
numerical  solution  will  be  based  on  the  assumption  of  local  reaction,  and  all  waves 
propagating  into  the  pores  will  be  assumed  to  propagate  normally  to  the  surface. 

Using  the  linearized  version  of  Attenborough’s  equations  (5.12)  to  (5.14)  with 
p'  eliminated,  and  assuming  one-dimensional  propagation  in  the  z  direction  with  z 
direction  particle  velocity  v,  the  time- dependent  equations  to  be  solved  numerically 
are 


c  dp’  dv  n 

ot  oz 

K  dv  *  ,  2  dp'  A 

+  +  Cq—z—  =  0 

uZ 


Po  n  at 


(9.1) 

(9.2) 


Recall  from  the  last  chapter  that  the  required  numerical  boundary  conditions  along 
the  ground  boundary  were  on  p'  and  v.  In  solving  the  air  and  porous  equations  to¬ 
gether,  the  porous  surface  will  lie  on  z  —  z^.  The  obvious  boundary  conditions 
for  this  case  are 


P  air  P  porous  (9.3) 

P sdr  ~  ^porous  (9.4) 

along  the  boundary  where  p'  ^  and  v ^  are  the  variables  in  Chapter  4  Equations 
(4.34)  and  (4.36)  and  p'^rous  and  ^porous  816  variables  in  (9.1)  and  (9.2).  Note 
that  although  Uporoug  is  an  average  particle  velocity  and  v^T  is  not,  that  (9.4)  is  still 


correct  because  they  are  both  quantities  defined  at  points  in  space. 
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To  solve  (9.1)  and  (9.2),  one  variant  of  the  second-order  in  space  and  time 
MacCormack  finite  difference  method  was  used,  giving 


n'ti  _  ~>n  P°  ^  i  „n\ 

-  ir^r  j+'  j 


fi  co 

—  „n _ 0 _ OL. 


„>n  + 1  _  1  i  Po  ^  ,,n  \ 

pJ  "T  pJ  +  pJ  7T3*  '  vJ~l) 


u/+1 


1  c02  A* 

-  ¥  +  *  -  7^'^  - ^  -  * 


(9.5) 


(9.6) 


(9.7) 


(9.8) 


This  finite  difference  scheme  must  follow  stricter  stability  guidelines  than  did  the 
equations  in  air,  due  to  the  presence  of  the  flow  resistivity  4>.  One  can  find  the  res¬ 
trictions  on  the  method  by  the  standard,  but  tedious,  procedure  of  performing  a 
Von  Neumann  stability  analysis.114  Let 


(9.9) 


then  one  can  determine  the  amplification  matrix  G  where 


An+1  An 

gn  +  l  ~  G(Atfkz)  gn 


(9.10) 


by  manipulating  Equations  (9.5)  to  (9.8).  The  result  for  low  frequencies,  kz  — ►C,  is 


1  0 

G<  Af,0)  -  Q  ^ 


(9.11) 


where 
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e  =  i- 


P0K 


n 


(9.12) 


The  eigenvalues  of  GK  A£,0)  are 


X 


1,2  = 


1 

2 


1  +  C±  (l+02-4e 


(9.13) 


from  the  solution  of  det(\I  —  GKAi.O))  =  0.  The  Von  Neumann  necessary  condi¬ 
tion  for  stability  states  that 


X, 


<  1  +0(At) 


(9.14) 


for  all  i.  Under  the  condition  -+0,  no  flow  resistance,  £  =  0  which  implies 


X, 


=  1,  and  stability  is  insured.  However,  for  0  one  can  see  from 


(9.12)  that  the  time  step  At,  relative  to 


P0K 


,  determines  whether  the  criterion 


(9.14)  is  met  It  was  found  that  choosing  a  At  too  large  in  the  porous  medium  nu¬ 
merical  solution  will,  in  fact,  make  the  numerical  solution  go  unstable. 

Since  at  low  frequencies  the  speed  of  sound  in  the  air  is  on  the  order  of  10 
times  faster  than  the  speed  of  propagation  in  a  porous  medium,  the  wavelengths  are 
on  the  order  of  10  times  shorter  in  the  porous  medium  than  in  the  air.  Because  of 
this  difference,  the  Ax  grid  spacing  in  the  pores  should  be  on  the  order  of  1/  10th  of 
the  size  of  the  grid  spacing  Ax  in  cir.  The  question  then  becomes  how  the  nu¬ 
merical  scheme  in  the  air  will  be  mated  with  the  numerical  scheme  in  the  pores, 
considering  this  difference  in  grid  spacings  in  the  air  and  pores. 


It  was  chosen  to  numerically  couple  the  pore  solution  to  the  air  solution  by 
refining  the  air  grid  in  the  z  direction  as  it  comes  into  contact  with  the  pore  grid. 
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Set  Ax^  =  lOAx^,,,  and  At^  =  lOAi^g.  Then,  referring  to  Figure  9.1,  refine 
the  air  between  j  =2  and  j  —  1  and  t={n+l)At  and  t=nAt  by  10.  Calculate  the  7=2, 
t=(n  + 1)^^,.  point  by  the  interior  scheme  using  Ax.^  and  Af^  from  points  7=1,  2, 
and  3  at  time  t=nAtair.  Then  interpolate  in  time  the  values  between  t—{n  +1)^^ 
and  n  At  a,  for  j— 2  to  obtain  the  values  by  the  time  increments  of  for  j=  2, 

indicated  by  the  x’s  in  Figure  9.1.  Next  calculate  the  values  of  air  at  the  spatial 
points  between  j=  2  and  7=1  and  the  values  of  the  porous  medium  at  all  spatial 
points  less  than  z  =  zmin,  using  the  refined  values  Ax^^  and  At p^,  until  the  7=1 
t={n+\)  At ^  values  are  calculated.  This  is  the  boundary  condition  point  needed  for 
the  next  iteration.  Thus,  air  values  between  7=1  and  7=2  and  the  pore  values  are 
solved  ten  times  as  often  as  the  air  values  in  the  rest  of  the  air  grid. 

In  1986  Marsha  Berger  analyzed  this  method  of  interpolating  in  time  and 
refining  in  time  and  space  thoroughly  for  stability.116  Her  analysis  is  for  a  model 
equation  solved  using  the  Lax-Wendroff  method,  very  similar  to  the  MacCormack 
techniques  used  in  this  thesis.  She  found  that  refining  in  both  time  and  space  by  in¬ 
teger  multiples  was  stable.  Note  that  if  a  nondissipative  difference  method  were 
used  for  the  computations  with  this  type  of  grid  refinement,  the  numerical  solutions 
may  not  be  stable. 

To  test  this  coupling  of  the  air  to  the  porous  ground,  computations  were  per¬ 
formed  to  estimate  the  actual  surface  impedance  of  the  porous  model.  It  is  possible 
to  determine  impedance  using  impedance  tube  techniques.116,117  Both  50  Hz  and  25 
Hz  linear  plane  wave  sinusoids  were  sent  propagating  in  the  — z  direction,  normally 
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< - air  pores - ) 


Figure  9.1.  Finite  difference  implementation  of  the  air-porous  interface.  All  the 
points  to  the  left  of  z  —  zmin  are  air  and  those  to  the  right  are  porous.  The  air 
points  between  7=2  and  7=1,  like  all  the  porous  points,  are  refined  in  z  and  t  by  10 
compared  to  the  rest  of  the  air  grid. 


Figure  9.2.  Impedance  measurement  maximum  pressure  over  time  as  a  function  of 
distance  from  the  air-porous  interface.  The  first  minimum  pmax  from  the  interface 
is  at  zn ,  and  the  first  maximum  is  at  zt . 
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onto  the  porous  medium  from  the  air.  A  standing  wave  pattern  results  from  the 
direct  wave  and  the  reflected  wave  superposing.  Figure  9.2  plots  the  maximum 
pressure  over  time,  pmax,  as  a  function  of  distance  from  the  air-porous  boundaiy  at 
z min.  The  particular  wave  envelope  shown  is  typical  of  the  type  seen  with  a  large 
standing  wave  ratio,  characteristic  of  a  high  magnitude  of  impedance.  To  calculate 
numerically  the  impedance,  the  three  numbers  needed  are  the  distance  zn  —  zmm, 
the  maximum  pmax  at  zx,  and  the  minimum  pmax  at  zn. 

The  measurements  made  were  difficult  and  crude  since  the  finite  difference 
grids  could  not  be  made  fine  enough  using  the  two-dimensional  program.  Finer 
grids  would  have  resolved  the  distance  of  the  minimum  pmax  to  the  air-pore  inter¬ 
face  more  precisely  but  would  have  used  a  tremendous  amount  of  computer  time. 
The  value  of  Ax  for  both  simulations  presented  here  was  0.0625  m,  and  the  peak 
pressure  of  the  sinusoid  was  140  dB,  making  the  wave  linear. 

For  25  Hz,  the  measured  value  of  the  impedance  was  Z/p0c0  = 
14.737  +  i33.146,  and  for  50  Hz,  the  value  was  Z  /p0c0  =  11.031  +  136.348.  These 
impedances  are  low  in  magnitude  compared  to  the  results  of  Attenborough’s  low 
frequency  interpretation  of  the  phenomenological  porous  model,  as  seen  in  Table 
5.1.  The  phases  are  significantly  different  also.  Unfortunately,  it  is  difficult  to  tell 
whether  the  difference  between  the  values  of  impedances  is  due  to  measurement  er¬ 
ror,  which  was  a  problem,  or  due  to  fundamental  differences  between 
Attenborough’s  predictions  and  the  numerical  results  for  the  model  porous  equa¬ 


tions. 
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However,  it  is  clear  that  the  values  found  by  the  numerical  impedance  tube 
technique  used  here  are  certainly  close  enough  so  that  one  can  infer  that  the  air- 
porous  boundary  is  working  correctly.  Real  impedance  tube  data  taken  from  out¬ 
door  measurements  have  quite  a  bit  of  scatter,  and  the  values  measured  numerically 
here  are  within  the  tolerances  seen  in  such  physical  experiments. 
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10.  VERIFICATION  of  the  numerical  solution 

This  chapter  demonstrates  the  validity  of  the  numerical  model  developed  in  the 
preceding  chapters.  Several  different  types  of  tests  were  performed.  First,  the 
results  of  the  numerical  solution  method  described  in  the  preceding  chapters  are 
compared  to  those  of  the  Pestorius  algorithm,  for  the  case  of  one-dimensional  pro¬ 
pagation  of  electric  spark  pulses  in  the  free  field.  This  example  will  show  that  the 
ratio  c0  Jg/Ax  must  be  about  0.25. 

Next,  spark  pulse  reflection  numerical  results  are  given  for  bom  normal  and 
oblique  incidence  from  a  hard  surface.  The  normal  reflection  simulation  is  com¬ 
pared  quantitatively  with  an  analytic  result,  while  the  oblique  reflection  numbers  are 
compared  qualitatively  with  the  effects  observed  for  explosions  on  a  larger  scale. 

In  addition,  the  influence  of  artificial  viscosity  on  the  numerical  solutions  is  ad¬ 
dressed,  for  both  the  cases  of  spark  pulse  and  blast  modeling.  Finally,  some  com¬ 
ments  are  made  on  the  value  of  /rB  used  for  the  blast  wave  runs  of  the  next 
chapter. 

1G.1  Electric  Spark  Pulses  in  the  Free  Field 

Spherical  spark  pulses  in  air  were  used  to  compare  one-dimensional  free  field 
propagation  using  the  algorithms  of  this  dissertation  and  the  Pestorius  algorithm. 
The  absorption  parameters  of  both  programs  corresponded  to  the  ANSI  standard  for 
attenuation  in  the  atmosphere.118  The  particular  version  of  the  Pestorius  algorithm 
used  was  that  of  Bass,  Ezell,  and  Raspet119  Although  this  program  adds  the  effects 
of  vibrational  relaxation  dissipation,  and  the  current  implementation  of  the  equa- 
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dons  of  Chapter  4  does  not,  the  frequencies  included  in  the  spark  pulse  were  high 
enough  to  not  be  affected  by  the  02  and  N2  relaxation  effects.  For  all  the  spark  cal¬ 
culations,  yug  was  set  to  zero,  commonly  known  as  the  Stokes  assumption.120 

The  values  /»0  =  1.21  kg/m3  and  c0  —  343  m/s  were  used  with  a  grid  322  x  322 
bounding  0.005  to  0.045  m  in  d  and  -0.02  to  0.02  m  in  z  where  Ax  =  1.25x  10~6  m 
m  the  finite  difference  program.  The  constant  Ai  was  picked  as  2/  =  100x  10~9  s, 
giving  c0  Ai/Ar  =  0.2744,  appropriate  for  the  algorithms  of  Chapter  7.  The  pulse 
was  started  with  the  initial  peak  pressure  of  ^)0  =  20,000  Pa  or  180  dB  referenced 
to  20  fu Pa  at  0.03  m  from  d  =  z  =  0.  The  values  of  t+,  t,  tx,  f2,  and  /3  were  given 
in  Chapter  6,  r  =  40x  10~6  s,  for  example.  The  pulse  was  propagated  for  400x  10-6 
s,  corresponding  to  4000 At. 

To  follow  the  pulse,  the  windowing  capability  of  the  program  was  employed  in 
the  d  direction.  Numerical  receivers  were  placed  at  z  =  0.0  and  d  =  0.03, 
0.057447,  0.08488,  0.11232,  and  0.13976  m  to  correspond  to  the  placement  of  the 
receivers  in  the  Pestorius  algorithm.  An  artificial  viscosity  coefficient  of  v  ~  0.08 
was  used,  and  the  simulation  took  18.56  min  on  a  Cray-2  supercomputer. 

The  solid  lines  in  Figures  10.1  to  10.5  show  the  results  of  this  run  while  the 
dashed  lines  provide  the  results  of  the  Pestorius  algorithm.  Figure  10.1  gives  the 
initial  conditions.  Note  that  the  initial  conditions  for  the  two  programs  are  not  ex¬ 
actly  the  same,  therefore  their  propagation  should  not  be  exactly  alike.  As  Figures 
10.2  to  10.5  show,  the  pulse  shape  becomes  closer  to  the  letter  N  as  it 


Figure  10.1.  A  free  field  spherical  spark  pulse  propagating  in  air.  The  solid  lines  are 
the  finite  difference  results,  and  the  dashed  lines  are  from  the  Pestorius  algorithm. 
Time  is  given  relative  to  a  wave  traveling  at  343  m !  s.  Initial  conditions  at  d  =  0.03 
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Figure  10.7.  Same  as  Figure  10.1  at  d  =  0.13976  m,  with  c0&/&c  =  0.5488  and 
no  artificial  viscosity. 
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propagaies.  This  is  the  well  kno  vn  nonlinear  phenomenon  of  propagation  speed  be¬ 
ing  a  function  of  amplitude.  The  positive  acoustic  pressures  travel  faster  than  the 
sound  speed  c., ,  while  the  negative  acoustic  pressures  crave!  slower  than  c0,  accord¬ 
ing  to  Equation  ‘2.10  of  Chapter  2.  Most  of  tne  nonlinear  steepening  occurs  eariy 
m  the  propagation  smee  the  distortion  is  a  function  o  '  amplitude  which  is  failing  off 
;  as  i/r. 

It  is  obvious  that  the  agreement  between  the  finite  difference  program  and  the 
Pestorius  algorithm  is  good.  The  pulses  lengthen  at  about  the  same  rate,  and  the 
amplitude  dependence  is  similar.  Thi ,  indicates  that  the  amounts  of  nonlinearity 
and  absorption  seen  in  tne  two  programs  were  nearly  the  same.  Notice,  however, 
that  the  finite  difference  program  takes  into  account  all  the  second-order  nonlinear 
effects,  while  the  Pestorius  algorithm  uses  only  the  P  plane  wave  nonlinear  effects. 
Thus,  the  spherical  pulse,  started  in  the  free  field,  should  not  propagate  exactly  as 
the  Pestorius  algorithm  dictates.  This  may  account  for  some  of  the  difference 
between  the  two  results.  The  negative  pressure  dip  before  the  steep  lise  in  pressure 
at  the  front  of  the  pulse  in  the  finite  difference  program  case  is  a  direct  result  of  the 
phase  dispersion  mentioned  in  Chapter  7. 

This  phase  dispersion  increases  if  the  artificial  viscosity  is  removed.  Figure  10.6 
shows  the  same  finite  difference  result  at  d  =  0.13976  m  «.v  does  Figure  10.5,  but 
with  v  =  0.  Clearly,  phase  dispersion  produces  undesirable  ripples.  As  previously 
stated,  this  is  due  to  the  nonlinear  effects  generatirg  frequencies  that  are  too  high  to 
be  well  resolved  by  the  spatial  grid. 
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With  i'  set  co  0,  one  can  also  see  the  effect  of  having  c0  2//2ir  too  high.  Figure 
10.7  gives  the  results  of  another  run  at  d  =  0.13976  m,  with  the  difference  that  2/ 
was  doubled  while  u  =  0.  For  this  value  of  A t ,  c0  Jtf  'Ax  Ls  0.549.  The  large  ripples 
o>  .ceding  the  front  of  the  spark  pulse  are  due  to  the  errors  of  the  time  discretization 
dominating  the  eno  -  ~f  the  spatial  discretization.  It  was  found  that  artificial  viscos¬ 
ity  does  little  to  help  give  a  clean  solution  when  cn  Ai/2u  is  too  high. 


10.2  Normal  Spark  Pulse  Reflection  from  a  Hard  Surface 


This  section  gives  the  results  of  the  program  for  finite  amplitude  electric  spark 
pulse  normal  reflection  from  an  infinite,  flat,  hard  surface  boundary.  An  analytic 
result  exists  fi  r  this  case  which  will  be  used  for  compazison. 

For  normal  plane  reflection  from  a  hard  surface,  Pfriem  has  derived  the  pres¬ 
sure  amplification  fcctor 


a 


A 

Po 


1 


1 


(10.1) 


where  p0  is  atmospheric  pressure,  p,  is  the  incident  total  pressure  (atmospheric  pres¬ 
sure  plus  acoustic  pressure),  and  =  27/(7  —  D  where  7  is  the  ratio  of  specific 
heats.121  Here  a  gives  the  ratio  of  the  acoustic  pressure  near  a  hard  surface  when  a 
wave  reflects  to  the  acoustic  pressure  with  the  surface  absent  as  the  wave  passes  by. 
The  expression  is  derived  under  the  assumption  of  an  ideal  gas,  and  is  easy  to  com¬ 
pute  since,  for  air,  \  =  7.  In  the  linear  limit  a  equals  2,  the  familiar  result  of  pres¬ 
sure  doubling. 
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By  making  two  runs  of  the  finite  difference  program,  one  with  a  hard  surface 
and  one  without,  the  ratio  of  pressures  was  found  and  a  predicted.  The  program 
was  run  in  the  free  field  for  the  initial  peak  pressures  of  154,  160,  168,  174,  and  180 
dB  on  the  domain  5.005  <  d  <  5.085  and  -0.02  <  z  <  0.02  m.  The  peak  of  the 
pulse  was  placed  5.03  meters  from  the  origin  initially  tc  make  the  propagation  nearly 
pmne.  Here,  all  the  computational  boundaries  used  the  absorbing  or  characteristic 
conditions  of  Chapter  8.  Then,  for  the  same  initial  conditions,  the  domain  5.005  < 
d  <  5.045  and  -0.02  <  z  <  0.02  meters  was  used  where  the  boundary  d  =  dmax  = 
5.045  meters  was  changed  to  be  hard.  This  hard  boundary  was  implemented  in 
finite  difference  notation  as 
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for  all  j  where  1  <i  <  I. 

Condition  ( 10.2)  merely  represents  the  continuity  of  density,  or  pressure  in  the 
linear  limit  Equation  ( 10.3)  implies  that  the  hard  wall  does  not  move  as  the  wave 
impinges  on  it  Equation  ( 10.4)  is  not  required  physically,  but  is  required  numeri¬ 
cally  due  to  the  Strang  method  of  splitting  in  the  numerical  technique  of  Chapter  7. 
Lastly,  condition  ( 10.5)  represents  the  continuity  of  T\  from  the  equations  of  state 


( 4.9)  and  ( 4.10). 


99 


Both  162  x  162  ( Ax  =  250x  10~6,  M  =  200x  10'9,  and  v  =  0.01)  and  322  x 
322  (Ax  =  125x  10~6,  A£  =  100x  1CT9,  and  v  =  0.08)  runs  were  used  in  the  predic¬ 
tions  of  a .  Figure  10.8  gives  the  Pfriem  result  as  a  solid  line  showing  a  as  a  func¬ 
tion  of  the  free  field  incident  peak  pressure  in  dB  referenced  to  20/i  Pa  The  filled 
circles  in  this  figure  are  the  results  for  the  162x  162  runs  of  the  program,  and  the 
filled  squares  are  for  the  322x  322  runs.  The  programs  measured  the  peak  pressures 
encountered  at  z  =  0  and  d  =  5.045  —  ( Ax/2)  m  when  the  hard  surface  was  present 
or  absent,  and  these  numbers  were  compared  to  obtain  a.  Notice  that  the  peak 
pressures  in  the  free  field  and  hard  surface  runs  do  not  occur  at  the  same  times. 
The  peaks  of  the  finite  amplitude  waves  travel  at  speeds  proportional  to  their  ampli¬ 
tudes,  which  are  different  for  the  two  runs. 

The  agreement  between  the  Pfriem  result  and  the  finite  amplitude  result  is  gen¬ 
erally  very  good.  This  is  despite  the  fact  that  the  Pfriem  result  assumes  lossless  pro¬ 
pagation,  and  the  finite  difference  calculation  accounts  for  classical  absorption 
effects.  The  322x  322  predictions  seem  to  deviate  from  the  Pfriem  result  as  the  in¬ 
cident  peak  pressure  approaches  180  dB.  It  is  difficult  to  tell  how  much  of  this  man¬ 
ifestation  is  due  to  the  phase  dispersion  of  the  numerical  scheme  and  how  much  is 
due  to  the  second-order  lossy  nonlinear  equations  from  Chapter  4  beginning  to 
break  down  at  the  high  reflected  pressures,  on  the  order  of  186  dB.  If  the  problem 
is  due  to  phase  dispersion,  one  should  be  able  to  obtain  better  agreement  by  adjust¬ 
ing  the  artificial  viscosity  and  by  using  a  finer  grid. 
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Figure  10.8.  Normal  reflection  amplification  factor  versus  free  field  incident  peak 
pressure  in  dB  referenced  to  20  n  Pa  for  sparks.  The  solid  line  is  the  Pfriem  result 
The  circles  and  squares  are  the  finite  difference  results  for  a  coarse  ( 162x  162)  and  a 
finer  ( 322x  322)  grid  run,  respectively. 
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10.3  Oblique  Spark  Pulse  Reflection  from  a  Hard  Surface 

Notice  that  the  verification  calculations  presented  thus  far,  free  field  spherical 
propagation  and  plane  wave  normal  reflection,  could  have  been  performed  with  a 
one-dimensional  program.  There  would  have  been  a  significant  decrease  in  CPU 
time  and  memory  usage  for  a  one-dimensional  code.  However,  a  two-dimensional 
nonlinear  effect  is  now  presented,  the  cumulative  interaction  of  a  spherical  electric 
spark  pulse  obliquely  incident  on  a  hard  surface,  which  cannot  be  computed  using  a 
one-dimensional  program.  There  is  no  analytic  solution  to  compare  with  the  finite 
difference  results  in  this  case,  so  we  will  examine  the  qualitative  agreement  between 
the  program  results  and  what  is  expected  from  physical  arguments  and  explosion 
data 

In  the  normal  incidence  plane  wave  case,  the  Pfriem  equation  and  the  numeri¬ 
cal  results  gave  a  ranging  from  2.00  to  about  2.15  depending  on  the  incident  ampli¬ 
tude.  For  normal  reflection,  the  incident  and  reflected  components  of  the  pulse  in¬ 
teract  over  the  pulse  transit  time.  For  oblique  incidence,  however,  the  incident  and 
reflected  components  of  the  pulse  can  interact  over  much  longer  periods  of  time. 
Therefore,  it  is  reasonable  to  expect  that  for  oblique  reflection  the  interacting  in¬ 
cident  and  reflected  components  of  the  pulse  could  produce  a  larger  or  smaller  than 
is  the  case  with  normal  reflection.  The  amplification  factor  a  would  be  smaller  if 
the  oblique  incident  and  reflected  waves  destructively  interfered  with  each  other 
more  than  in  the  normal  incidence  case,  and  a  would  be  larger  if  they  constructively 
interfered.  It  also  stands  to  reason  that  the  interference  producing  deviations  in  a 
would  be  a  function  of  the  incident  angle. 
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Researchers  have  observed  this  phenomenon  experimentaily  and  numerically 

for  spherical  explosion  blasts  incident  on  natural  ground  surfaces.122,12,1,124  For  ex¬ 
plosions,  a  is  the  normal  reflection  result  for  angles  near  normal,  increases  to  a 
peak,  and  then  actually  decreases  below  2  for  angles  close  to  grazing  incidence.  The 
specific  curve  of  a  versus  angle  depends  on  the  incident  peak  pressure. 

Calculations  were  made  to  show  that  spark  pulses  have  a  versus  angle  curves 
similar  to  those  for  explosions.  The  program  was  run  for  a  spherical  spark  pulse 
above  a  hard  surface  atz  =  zmin  using  the  finite  difference  relations 
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for  all  i  where  1  <  t  <  /.  These  equations  are  similar  to  ( 10.2)  to  ( 10.5). 

The  spark  pulses  were  always  centered  around  the  origin  of  the  coordinate  sys¬ 
tem  initially,  and  the  hard  surface  was  located  atz^  =  -0.0075,  -0.01125,  -0.015, 
-0.01875,  -0.0225,  -0.025,  -0.03,  and  -0.035  m.  Corresponding  free  field  runs  on  a 
larger  domain  with  the  hard  surface  absent  were  also  made.  Numerical  receivers 
were  placed  along  the  ground  surface  or  its  equivalent  to  obtain  numerical  predic¬ 
tions.  From  one  pair  of  free  field  and  hard  surface  runs  it  was  possible  to  obtain 
one  point  to  put  on  each  plot  of  a  versus  incident  angle  for  each  constant  incident 
pressure.  Thus,  the  8  pairs  of  runs  produced  up  to  8  points  on  each  plot  of  a 
versus  incident  angle  for  fixed  incident  pressure. 
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Figure  10.9.  Oblique  reflection  amplification  factor  versus  angle  in  degrees  for  a 
fixed  free  field  incident  peak  pressure  for  sparks.  When  the  angle  is  0  degrees,  the 
wavefront  is  locally  normally  incident,  and  when  the  angle  is  90  degrees,  the  wave- 
front  is  locally  grazing.  The  free  field  incident  peak  pressure  is  15.0  kPa,  177.5  dB. 


Figure  10.10.  Same  as  Figure  10.9  except  the  peak  pressure  is  10.0  kPa,  174  dB. 
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Figure  10.11.  Same  as  Figure  10.9  except  the  peak  pressure  is  5.0  kPa,  168  dB. 


Figure  10.12.  Same  as  Figure  10.9  except  the  peak  pressure  is  2.5  kPa,  162  dB. 
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Figure  10.13.  A  composite  of  Figures  10.9  to  10.12.  Oblique  reflection 
amplification  factors  versus  angle  for  the  four  free  field  incident  peak  pressures  for 
sparks  of  15.0,  10.0,  6.0,  and  2.6  kPa 
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These  simulations  used  Ax  =  250x  10~°  m,  At  =  200x  10-9  s,  and  u  =  0.01 
with  the  initial  peak  0.03  m  from  the  origin  at  a  pressure  of  180  dB.  Each  of  the  16 
runs  needed  between  2.34  and  5.76  min  of  CPU  time  on  a  Cray-2. 

The  individual  data  points  in  Figures  10.9  to  10.12  present  the  results  of  these 
runs,  corresponding  to  peak  incident  pressures  of  15G00  Pa  (177.5  dB),  10000  Pa 
( 174  dB),  5000  Pa  ( 168  dB),  and  2500  Pa  (162  dB).  Figure  10.13  joins  each  of  the 
numerical  data  points  in  Figures  10.9  to  10.12  and  plots  these  curves  for  comparis¬ 
on.  The  results  confirm  the  assertion  that  a  should  be  a  function  of  incident  angle 
for  each  constant  incident  pressure.  As  can  be  seen  for  the  5000  Pa  numbers  (Fig¬ 
ure  10.11),  a  is  near  the  Pfriem  normal  incidence  result  for  most  small  angtes  of  in¬ 
cidence.  The  amplification  factor  a  increases  to  2.41  as  the  incident  angle  increases, 
and  then  it  decreases  below  2  to  1.78  or  less  as  the  angle  becomes  closer  to  grazing. 

These  curves  are  similar  to  those  seen  from  explosions  on  a  larger  scale  as 
given  in  the  Glasstone,  Kinney  and  Graham,  and  Heaps  et  al.  references.  This 
qualitative  agreement  between  the  spark  and  blast  oblique  incidence  a  versus  in¬ 
cident  angle  curves  is  intuitively  satisfying,  but  some  laboratory  work  should  be  per¬ 
formed  in  the  future  to  experimentally  verify  these  spark  oblique  curves.  No  refer¬ 
ences  to  this  type  of  experimental  data  have  been  found  in  the  literature. 

10.4  Artificial  Viscosity  and  Bulk  Viscosity  Values  for  Blasts 

Recall  that  u  =  0.08  was  used  for  the  free  field  spark  simulations  to  produce 
Figures  10.1  to  10.5.  This  value  was  found  by  trial  and  error,  increasing  and  de¬ 
creasing  until  the  best  result  was  obtained.  The  goal  was  to  have  the  minimum  u 
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which  would  still  reduce  the  ripples  produced  by  the  phase  dispersion. 

Numerous  calculations  were  made  with  idealized  blast  pulses  of  approximate 
length  t~  40x  10-3  s  and  various  peak  pressures.  Computations  were  performed  re¬ 
peatedly  to  find  optimum  values  of  u.  It  was  deemed  that  a  value  of  u  was  good 
when  it  was  the  minimum  value  such  that  the  drop  from  peak  pressure  to  the  zero 
crossing  was  monotonic,  when  no  two  values  of  time  between  the  peak  and  the  zero 
crossing  were  associated  with  the  same  value  of  pressure.  The  values  that  were 
found  for  =  200x  10-6  s  and  c0  £t/&x  =  0.2744  are  given  in  Table  10.1,  and  the 
values  for  -  100x  10“6  with  the  same  c0  At/Ax  in  Table  10.2. 

These  tables  demonstrate  that  the  amount  of  artificial  viscosity  needed  general¬ 
ly  increases  as  the  peak  amplitudes  increase.  This  makes  sense,  since  more  high 
frequencies  that  the  grid  cannot  well  resolve  are  generated  as  the  amplitude  in¬ 
creases.  More  nonlinearity  means  more  high  frequencies,  so  a  higher  u  is  needed 


Table  10.1.  Coarse  grid  minimum  artificial  viscosities  for  blasts. 
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for  higher  amplitudes.  The  fact  that  a  higher  u  is  needed  for  the  finer  grid  also 
makes  sense  from  Equation  (7.24)  of  Chapter  7.  That  relationship,  however,  is  not 
exactly  followed,  so  other  effects  not  directly  considered  here  may  be  involved. 

Finally,  the  values  of  /.i g  used  in  the  blast  runs  above  should  be  discussed. 
The  lossy  second-order  nonlinear  model  of  Chapter  4  -  arrently  does  not  include  the 
effects  of  molecular  relaxation  dissipation.  Relaxation  is  not  important  for  spark 
pulses  due  to  their  high  frequency  content,  so  Hg  was  assumed  negligible  and  set  to 
zero  for  sparks.  For  the  lower  frequencies  involved  with  blasts,  however,  one  can 
increase  fig  to  account  for  some  of  the  energy  lost  by  relaxation  processes. 

It  is  possible  to  derive  an  expression  for  the  absorption  due  to  a  relaxation  pro¬ 
cess  in  the  low  frequency  limit,  as  an  increment  to  Hg,  iVa-  The  expression  is 


-Va  - 


rc 


(  Ct:\) 


(10.10) 


where  rt  is  the  ith  relaxation  time  and  (a^X)m  is  the  maximum  absorption  per 
wavelength  associated  with  the  ith  relaxation  process.126 

If  one  looks  at  a  chart  of  absorption  coefficient  a  versus  f  for  air  such  as  Fig¬ 
ure  10-13  in  Pierce’s  book,  it  is  obvious  that  two  relaxation  processes  in  air,  02  and 
N*  primarily  control  absorption  at  low  frequencies.  Unfortunately,  using  the  low 
frequency  limit  for  both  these  processes  would  make  the  absorption  too  high  for 
frequencies  in  the  range  100  to  1000  Hz.  Using  too  high  an  absorption  would 
severely  affect  the  peak  sound  pressure  levels  observed  at  long  ranges. 
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The  finite  difference  numerical  method  described  in  Chapter  7  already  has 
some  intrinsic  dissipation,  to  which  an  artificial  fourth-order  viscosity  was  added. 
To  keep  the  peak  sound  pressure  levels  as  accurate  as  possible,  the  minimum  addi¬ 
tional  bulk  viscosity  was  sought.  Such  peak  levels  are  dominated  by  the  high  fre¬ 
quency  absorption.  If  the  low  frequency  limit  of  the  02  relaxation  process  alone  is 
used,  the  absorption  for  the  high  frequencies  in  the  blast  pulse  turns  out  to  be  about 
right,  and  this  technique  was  employed.  Oxygen  has  a  relaxation  frequency  at  ft  = 
12,500  Hz. 

For  the  02  bulk  relaxation  effect  (atX)m  is  0.0011  and  rt  =  1/(2 = 
12.732x  10-6  and  is  1.27x  10-3  using  c0  =  343  and  p0  =  1.21.  Although  pB 

is  usually  assumed  zero,  experimental  evidence  has  shown  it  to  have  the  value  // B 
=  0.6/j  for  air.  To  this  value  was  added,  giving  pB  =  0.6yu  +  \iB  = 


1.281x  10  3,  which  is  about  70  p. 
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11.  FINITE  \MPLITUDE  BLAST  PREDICTIONS 

This  section  chc  bindings  of  this  dissertation  on  the  propagation  of 

unite  amplitude  blast  sounds  outdoors.  Throughout,  the  emphases  is  on  the  oblique 
renecnon  of  bum.  sounds  from  hard  surfaces  and  from  -he  porous  medium  model  of 
„r,a  sun;;.  in  Chapters  5  and  S.  7  particular  ci*arge  size  and 

me  try  three  runs  were  made  Tne  first  run  was  in  the  free  field  with  no 
reflective  surface  present.  The  second  run  introduced  the  hard  surface  ground, 
while  the  third  used  the  porous  medium  ground.  This  chapter  will  hrst  give  the 
results  for  a  wide  variety  of  angles  of  Incidence,  corresponding  to  one  charge  at  vari¬ 
ous  heignts.  Then,  to  investigate  propagation  near  the  ground,  the  geometry  will  be 
held  fixed  for  differently  sized  charges. 

11.1  A  Fixed  Charge  Weight  at  Various  Heights. 

Numerous  runs  at  various  heights  were  made  that  would  have  had  a  peak  free 
held  pressure  of  180  dB  at  30  m  from  the  source.  Looking  at  Table  6.1,  a  180  dB 
pulse  at  30  m  is  183.89  ms  long  in  time,  or  equivalently,  approximately  63  m  in 
space.  This  spatial  dimension  implies  that  for  the  initial  conditions  of  this  disserta¬ 
tion  the  wave  cannot  be  within  63  m  of  the  origin.  For  the  runs  made  here,  the  ini¬ 
tial  condition  used  was  a  spherical  183.89  ms  pulse  that  had  a  peak  of  8,571.4  Pa  70 
m  from  the  origin.  This  new’  peak  makes  the  gross  linear  assumption  that  between 
30  m  and  /0  m,  the  wave  fell  off  as  1/r  and  did  not  distort.  These  calculations  were 
made  with  different  domains,  depending  on  the  heights  of  the  sources.  For  each 
run  the  artificial  viscosity  used  was  u  =  5x  10~3  with  a  coarse  grid  of  .At  =  200x  10“6 
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s  and  Ax  =  0.250  m.  The  three  programs  each  took  approximately  19  to  30  min  on 
a  Cray-2  for  each  height  investigated. 

Graphs  were  developed  giving  the  amplification  factor  a  as  a  function  of  in¬ 
cident  angle  for  a  fixed  incident  pressure.  As  in  the  last  chapter,  a  is  the  ratio  of 
the  acoustic  pressure  received  near  the  reflection  surface  to  the  pressure  at  the  same 
receiver  with  the  reflection  surface  removed.  The  source  heights  used  to  obtain  the 
numerical  predictions  were  1,  11,  21,  31,  41,  51,  and  61  m.  Figures  11.1,  11.2, 
11.3,  and  11.4  give  the  resulting  hard  surface  amplification  factor  versus  angle  plots 
for  the  fixed  incident  free  field  peak  pressures  of  9000,  5000,  2500,  and  1250  Pa, 
respectively.  Figure  11.5  presents  composite  graphs  when  all  the  discrete  data  points 
from  Figures  11.1  to  11.4  are  joined  by  lines. 

These  hard  surface  amplification  factor  versus  incident  angle  plots  show  two 
mqjor  results.  First,  the  amplification  factors  are  up  to  2.65  for  the  blast  runs, 
whereas  they  were  up  to  only  2.4  for  the  hard  oblique  runs  of  electric  spark 
reflections  in  the  last  chapter.  Secondly,  the  hard  surface  curves  of  a  versus  in¬ 
cident  angle  are  also  more  peaked  for  the  blasts  in  Figure  11.5  than  for  the  sparks  in 
Figure  10.13.  This  difference  in  hard  oblique  incidence  amplification  factors 
between  sparks  and  blasts  indicates  that  these  reflections  do  not  scale  simply. 

Similarly,  the  oblique  porous  run  results  for  amplification  factors  appear  in  Fig¬ 
ures  1 1.6  to  11.10.  In  general,  these  results  are  more  peaked  and  the  peaks  are 
more  toward  grazing  incidence  compared  to  the  hard  cases.  It  is  unwise  to  directly 
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Figure  11.1.  Oblique  hard  reflection  amplification  factor  versus  angle  in  degrees  for 
a  fixed  free  field  incident  peak  pressure  for  blasts.  When  the  angle  is  0  degrees,  the 
wavefront  is  locally  normally  incident,  and  when  the  angle  is  90  degrees,  the  wave- 
front  is  locally  grazing.  The  free  field  incident  peak  pressure  is  9.0  kPa,  173  dB. 
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Figure  112.  Same  as  Figure  11.1  except  the  peak  pressure  is  5.0  kPa,  168  dB. 


Figure  11.3.  Same  as  Figure  11.1  except  the  peak  pressure  is  2.5  kPa,  162  dB. 
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Figure  11.4.  Same  as  Figure  11.1  except  the  peak  pressure  is  1.25  kPa,  156  dB. 
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Figure  11.5.  A  composite  of  Figures  11.1  to  11.4.  Oblique  hard  reflection 
amplification  factors  versus  angle  for  the  four  free  field  incident  peak  pressures  for 
blasts  of  10.0,  5.0,  2.5,  and  1.25  kPa 
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Figure  11.6.  Oblique  porous  reflection  amplification  factor  versus  angle  in  degrees 
for  a  fixed  free  field  incident  peak  pressure  for  blasts.  When  the  angle  is  0  degrees, 
the  wavefront  is  locally  normally  incident,  and  when  the  angle  is  90  degrees,  the 
wavefront  is  locally  grazing.  The  free  field  incident  peak  pressure  is  9.0  kPa,  173 


Figure  11.7.  Same  as  Figure  11. o  except  the  peak  pressure  is  5.0  kPa,  168  dB. 
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Figure  11.8.  Same  as  Figure  11.6  except  the  peak  pressure  is  2.5  kPa,  162  dB. 
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3.5 


Figure  11.9.  Same  as  Figure  11.6  except  the  peak  pressure  is  1.25  kPa,  156  dB. 


Figure  11.10.  A  composite  of  Figures  11.6  to  11.9.  Oblique  porous  reflection 
amplification  factors  versus  angle  for  the  four  free  field  incident  peak  pressures  for 
blasts  of  10.0,  5.0,  2.5,  and  1.25  kPa 
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compare  the.  actual  amplification  factors  between  the  hard  and  porous  cases,  howev¬ 
er,  since  the  initial  conditions  for  the  two  cases  are  approximate. 

The  a  versus  angle  curves  developed  here  are  similar  to  the  strong  shock  calcu¬ 
lation  results  of  Heaps  et  al.  and  the  curves  in  Glasstone’ s  work,  The  Effect  of  Nu¬ 
clear  Weapons.  These  authors  do  not  use  the  notation  used  in  this  thesis,  however. 
For  comparison  Heaps  et  al.  call  a  the  reflected  overpressure/P,  while  Glasstone 
merely  refers  to  it  as  the  reflected  overpressure  ratio  Pr/Ps  where  Ps  is  given  in  at¬ 
mospheres.  The  similarity  of  the  curves  presented  in  this  thesis  to  the  curves  of  the 
Heaps  et  al.  and  Glasstone  references  shows  that  the  numerical  method  of  this 
thesis  is  accurately  predicting  the  oblique  incidence  of  blast  waves. 

11.2  Different  Charge  Sizes  at  Near  Grazing  Incidence. 

Many  rims  were  made  with  source  heights  of  1  m  with  various  initial  peak  pres¬ 
sures  to  see  the  effect  of  a  varying  charge  size  for  a  fixed  geometry.  Runs  with  ini¬ 
tial  peak  pressures  corresponding  to  180,  174,  168,  162,  156,  and  150  dB  at  30  m 
were  carried  out  The  180  dB  and  174  dB  cases  were  performed  using  initial  peaks 
of  8,571.4  Pa  at  70  m  and  7,500  Pa  at  40  m,  respectively,  since  these  larger  blasts 
could  not  be  started  at  30  m.  The  pulse  initial  durations  and  corresponding  charge 
sizes  of  C-4  plastic  explosive  are  listed  in  Table  5.1. 

Receivers  at  30,  45,  60,  90,  120,  180,  240,  360,  480,  720,  and  960  m  from  the 
source  at  heights  of  0,  1,  2,  and  5  m  recorded  the  waves  as  they  propagated,  as  well 
as  all  maximum  pressures  encountered  near  the  ground.  These  free  field,  hard  sur¬ 
face,  and  porous  surface  runs  were  made  using  the  artificial  viscosities  of  Table 
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10.2,  or  interpolated  values  thereof,  with  =  100x  10-6  s  and  c0^t/^x  =  0.2744. 
The  150  dB,  156  dB,  and  162  dB  peak  runs  were  propagated  for  0.4  or  0.2  s.  The 
168  dB  run  modeled  propagation  for  0.8  s,  the  174  dB  run  for  1.6  s,  and  the  180  dB 
run  for  3.2  s.  The  shorter  programs  took  about  20  min  on  a  Cray-2,  while  the  larg¬ 
est  programs  took  5  h  or  more  of  Cray-2  CPU  time  apiece.  The  180  dB  runs  were 
particularly  expensive  computationally  because  of  the  laige  physical  size  the  blast 
pulse  occupied,  and  their  computation  would  not  be  possible  using  something  other 
than  a  state  of  the  art  supercomputer  ( as  of  1990). 

The  resulting  150,  162,  and  174  dB  peak  SPL  versus  range  plots  for  propagation 
near  the  hard  surface  are  presented  in  Figures  11.11  to  11.13,  respectively.  In  these 
figures  the  0  m  receiver  height  peaks  are  given  by  the  solid  lines,  the  1  m  receivers 
by  the  short  dashed  lines,  and  the  5  m  receivers  by  the  long  dashed  lines.  The  156, 
168,  and  180  dB  results  were  similar.  The  curves  show  that  although  the  peak 
sound  pressure  levels  are  initially  different  at  the  different  heights,  this  goes  away 
quickly  and  there  is  no  difference  in  level  as  a  function  of  height  further  out  The 
peak  levels  all  fall  off  near  the  rate  of  r~12,  7.23  dB  per  doubling  of  distance,  which 
is  given  on  these  figures  as  a  solid  line  denoted  by  “ref.”  This  r-1 2  reference  line  is 
started  at  170  dB  at  20  m  and  extends  to  640  m. 

This  decay  rate  of  approximately  r-1-2  agrees  with  the  empirical  relation  given 
by  Reed  for  weak  shocks  in  the  far  field  of  a  strong  blast126  The  finite  difference 
simulation,  therefore,  is  making  predictions  consistent  with  the  weak  shock  region 
of  a  strong  shock  calculation. 
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Figures  .  11.14  to  11.19  show  the  hard  ground  numerical  receiver  pressure 
versus  time  curves  at  120  m  range  for  a  150  dB  peak  free  field  initial  pressure  at  a  0 
m  height,  150  dB  peak  at  5  m,  162  dB  peak  at  0  m,  162  dB  peak  at  5  m,  174  dB 
peak  at  0  m,  and  174  dB  peak  at  a  5  m  height,  respectively.  The  150  dB  plots  show 
much  dispersion  and  absorption,  while  the  162  dB  plo cs  show  less,  and  the  174  dB 
even  less  than  that  As  is  expected  the  “shock”  is  more  pronounced  at  the  higher 
amplitudes.  Little  difference  is  seen  in  the  waveforms  between  when  the  numerical 
receiver  is  near  ground  level  and  when  the  receiver  is  at  a  5  m  height 

The  porous  simulation  peak  SPL  versus  range  plots  for  150,  162,  and  174  dB 
are  given  in  Figures  11.20  to  11.22.  Again,  the  results  for  156,  168,  and  180  dB 
peak  sound  pressure  levels  were  similar.  For  this  porous  case  the  plots  show  a 
marked  difference  in  peak  levels  as  a  function  of  height,  compared  to  the  hard 
simulations  which  showed  little  difference.  The  peak  sound  pressure  levels  were  al¬ 
ways  the  greatest  near  the  porous  surface.  This  effect  may  seem  counterintuitive, 
but  is  due  to  the  near  grazing  geometry  used  here  and  the  particular  impedance  of 
the  porous  surface,  i.e.,  a  high  magnitude  impedance  with  phase  angle  of  7r/4. 

For  the  0  m  receiver  height  the  porous  peak  levels  decayed  at  about  r-1 2,  ex¬ 
cept  for  the  150  dB  initial  peak  run  which  fell  off  faster.  This  rate  of  r-1-2  is  what 
was  seen  for  the  hard  surface  runs.  However,  for  increasing  initial  peak  pressure, 
the  sound  pressure  levels  at  the  5  m  height  decayed  at  a  rate  less  than  r-1 2  and 
those  at  the  1  m  height  even  less  than  that.  Comparing  Figure  11.13  1 174  dB,  hard) 
to  Figure  11.22  ( 174  dB,  porous)  after  the  simulation  settles  down,  it  is  clear  that 
the  slopes  of  the  curves  for  the  0  m  height  cases  are  nearly  identical. 
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Figure  11.14.  Pressure  versus  time  at  the  120  m  range  and  the  0  m  receiver  height 
for  a  150  dB  initial  peak  blast  propagating  over  the  hard  surface. 
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Figure  11.17.  Pressure  versus  time  at  the  120  m  range  and  the  5  m  receiver  height 
for  a  162  dB  initial  peak  blast  propagating  over  the  hard  surface. 
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Here  one  also  sees  that  the  1  m  height  short  dashed  line  is  falling  off  significantly 
less  fast  than  the  solid  0  m  height  line,  and  the  5  m  height  long  dashed  line  slightly 
less  fast 

These  results  indicate  that  the  peak  amplitudes  are  decaying  less  fast  for  in¬ 
creased  charge  sizes  for  receivers  not  on  the  porous  ground  surface.  Larger  charge 
size  blasts,  hence,  seem  to  undergo  less  attenuation  due  to  ground  impedance  than 
do  smaller  charge  size  blasts.  Larger  blasts  have  higher  initial  peak  levels  to  in¬ 
crease  the  nonlinear  effects,  but  also  are  longer  and  have  an  increased  low  frequen¬ 
cy  content  Since  both  frequency  content  and  peak  amplitudes  are  being  changed 
for  different  size  charges,  the  decrease  in  “excess  attenuation”  for  larger  charges 
clearly  cannot  be  attributable  solely  to  the  increased  finite  amplitude  effects. 

The  pressure  versus  time  plots  for  the  porous  runs  are  shown  in  Figures  11.23 
to  11.28.  Here  the  150  dB  initial  free  field  peak  plots  at  heights  0  m  and  5  m.  Fig¬ 
ures  11.23  and  11.24,  respectively,  show  that  the  amplitude  magnitudes  of  both  the 
peak  and  rarefaction  are  significantly  larger  near  the  ground.  This  is  quite  different 
than  for  the  hard  surface  case  where  virtually  no  amplitude  difference  was  detected 
as  a  function  of  height  For  the  larger  initial  peak  amplitude  cases  in  Figures  11.25 
to  11.28,  however,  only  the  peak  of  the  blast  waves  showed  a  significant  increase  in 
pressure  near  the  ground.  Here  the  “trough”  shapes  and  magnitudes  of  the  pulses 
seem  to  be  unchanged  as  a  function  of  height  for  both  162  and  174  dB. 

These  findings  agree  with  time  harmonic  linear  theory  for  the  particular 
geometry  and  complex  impedance  used  here.  Sngle  frequency  waves  were  run  in  a 
homogeneous  atmosphere  sound  propagation  program  based  on  the  work  of 


160.0 


120  0 


80  00 


40  00 


0  000 


-40  00 


-80  00 


-120  0 


-160  0 


8  2 


66  .2640  2714  2788  2862  2936 


t  /  S 


Figure  11 
for  a  160 


the  120  m  range  and  the  0  m  receiver  height 


-  ?500  H - 1 - 1 - 1 - 1 - 1 - 1 - 1 - 1 

2289  2417  2544  .2672  2800  2928  3055  3183  3311 

t  ,  S 


Figure  11.25.  Pressure  versus  time  at  the  120  m  range  and  the  0  m  receiver  height 
for  a  162  dB  initial  oeak  blast  DroDa£ralin£r  over  the  norous  surface 
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Figure  11.26.  Pressure  versus  time  at  the  120  m  range  and  the  5  m  receiver  height 
for  a  162  dB  initial  peak  blast  propagating  over  the  porous  surface. 
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Nobile  and  Hayek.40’127  In  the  program  the  values  of  impedance  from  Table  5.1 
were  used,  as  a  function  of  frequency.  This  Nobile  and  Hayek  calculation  indicated 
that  there  was  little  difference  in  levels  as  a  function  of  height  for  low  frequency 
waves,  but  a  significant  difference  as  a  function  of  height  for  high  frequencies. 

In  Figures  11.23  to  11.23  the  150  dB  initial  peak  pulse  contains  more  high  fre¬ 
quency  energy  relative  to  higher  amplitude  blasts  which  allows  both  the  peak  and 
trough  of  the  pulse  to  be  greater  in  magnitude  at  the  0  m  height  than  at  the  5  m 
height  For  the  larger  amplitude  waves,  however,  the  high  frequencies  manifest 
themselves  near  the  “shock”  fronts  of  the  pulses  and  it  is  here  that  the  amplitude 
difference  is  seen  between  the  two  heights.  Another  contributor  to  the  peaks  being 
increased  is  certainly  finite  difference  phase  dispersion.  This  numerical  effect  would 
enhance  the  peaks  also,  particularly  for  the  highest  runs  made  at  174  and  180  dB  in¬ 
itial  free  field  pressures. 

Pom  us  and  hard  ground  runs  were  also  made  with  the  source  heights  of  3  m 
instead  of  1  m.  These  values  will  not  be  reported  in  detail  here  for  brevity's  sake, 
since  the  values  were  close  to  those  for  the  1  m  height  blasts.  No  difference  was 
seen  in  the  far  field  for  either  the  hard  or  porous  cases,  and  near  field  differences 
were  only  1  to  2  dB. 
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12.  CONCLUSIONS  AND  FUTURE  WORK 

This  dissertation  has  described  a  finite  difference  numerical  model  for  the 
propagation  of  finite  amplitude  acoustic  blast  waves  in  a  homogeneous  lossy 
atmosphere.  The  nonlinear  air  propagation  equations  included  classical  dissipative 
effects  and  may  easily  be  extended  to  include  molecular  relaxation  effects.  The 
model  also  accounts  for  both  hard  and  simple  porous  media  ground  surfaces.  A 
time-dependent  locally  reacting  formulation  of  Morse  and  Ingard’s  phenomenologi¬ 
cal  model  of  a  porous  medium  has  been  developed  to  represent  the  porous  ground. 
The  propagation  of  finite  amplitude  pulses  such  as  electric  spark  and  blast  waves  has 
served  to  validate  the  computer  algorithms. 

The  model  agrees  with  the  results  of  free  field  propagation  of  sparks  using  the 
Pestorius  algorithm.  It  further  agrees  with  the  analytic  predictions  of  the  pressure 
amplification  factor  a  for  the  case  of  normal  plane  reflection.  The  blast  computa¬ 
tions  of  a  versus  angle  for  various  fixed  incident  peak  pressures  also  give  similar 
results  to  those  models  cum.  ltly  in  use  for  strong  shock  waves.  The  finite 
difference  program  computes  slightly  different  a  versus  angle  curves  when  the 
bounding  surface  is  changed  from  hard  to  porous. 

For  different  charge  sizes  detonated  near  the  hard  ground  surface,  the  simulat¬ 
ed  peak  sound  pressure  levels  decayed  very  closely  to  the  power  law  r-1,2  which  has 
been  observed  in  physical  experiments.  It  was  found  that  the  larger  the  charge  size, 
the  less  the  peak  levels  above  the  ground  are  affected  by  finite  ground  impedance 
excess  attenuation.  This  is  an  important  result,  not  predicted  by  linear  theory  pro¬ 
pagation  codes. 
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The  interaction  between  finite  ground  impedance  and  finite  amplitude  effects 
seems  to  be  significant,  but  not  straightforward.  The  assumption  of  a  first-order  ad¬ 
ditive  combination  of  ground  impedance  and  nonlinear  effects  should  not  predict,  in 
general,  the  correct  decay  rates  of  peak  sound  pressure  levels  with  distance.  Such 
an  additive  interaction  between  finite  amplitude  sound  and  ground  impedance  was 
proposed  by  Raspet,  Bass,  and  Ezell,  for  example.128 

The  computer  program  may  be  easily  extended  to  increase  its  accuracy  and  use¬ 
fulness.  The  first  obvious  extension  to  the  code  is  to  add  molecular  relaxation 
effects  as  described  in  Chapter  4.  Using  such  relaxation  effects  in  a  one-dimensional 
code  would  allow  the  user  to  calculate  blast  wave  rise  times.  This  one-dimensional 
program  would  be  fast  compared  with  the  two-dimensional  version  and  would  run 
on  small  mainframe  computers  instead  of  on  a  supercomputer.  Another  useful  ex¬ 
tension  would  be  to  permit  p0  and  c0  to  vary  with  space.  This  addition  would  allow 
the  modeling  of  nonlinear  propagation  in  a  refractive  atmosphere. 

It  would  also  be  desirable  to  use  a  1-D  program  to  further  test  the  impedance 
of  the  air-ground  interface  presented  in  this  thesis,  over  a  wide  variety  of  frequen¬ 
cies.  Recall  that  computation  expense  made  surface  impedance  measurements 
difficult  for  the  2-D  program.  Since  the  amplitudes  of  the  waves  in  the  air  are 
finite,  it  also  makes  sense  that  a  better  porous  model  would  account  for  nonlinear 
effects.  One  might  use  the  nonlinear  flow  resistance  as  developed  by  Blackstock  et 
al  for  this  modification.  Another  avenue  to  pursue  to  give  a  more  realistic  porous 
medium  model  would  be  to  make  the  porosity  of  the  medium  vary  with  depth. 
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The  newly  developed  numerical  techniques  called  ENO  for  Essentially  Non  Os¬ 
cillatory  and  TVD  for  Totally  Variational  Diminishing  might  be  better  numerical 
techniques  for  solving  the  equations  of  this  research.129,130  The  ENO  methods  have 
only  recently  been  used  successfully  in  more  than  one  dimension,  and  questions 
remain  on  how  TVD  might  work  in  multidimensions.  For  a  1-D  code  these 
methods  should  eliminate  almost  all  of  the  phase  dispersion  seen  with  the  MacCor- 
mack  schemes. 

Another  area  that  merits  further  research  is  the  development  of  absorbing 
boundary  conditions.  It  would  be  useful  if  an  absorbing  boundary  not  only  exactly 
absorbed  energy  from  one  isolated  source,  but  also  from  a  source  and  its  images. 
The  recent  work  of  Gustaffson  provides  some  insight  into  possible  implementations. 
In  addition,  this  dissertation  contains  linear  lossless  boundary  conditions  which  were 
applied  to  a  system  of  lossy  nonlinear  equations.  One  should  derive  the  absorbing 
boundaries  for  the  lossy  nonlinear  system  itself  and  provide  a  mathematical  stability 
proof  for  the  resulting  conditions. 

Finally,  this  thesis  has  centered  solely  on  the  propagation  of  sound  in  air  out¬ 
doors.  The  equations  developed,  however,  are  applicable  to  other  fluid  acoustic 
media  also,  such  as  water.  There  is  great  promise  in  using  this  type  of  numerical 
solution  in  biomedical  ultrasonics,  for  example,  lithotripsy,  diagnostic  ultrasound, 
and  hyperthermia  modeling.  Applying  the  finite  difference  methods  in  the  time 
domain  to  underwater  acoustics,  such  as  in  simulating  parametric  arrays,  should  be 
possible  also. 
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APPENDIX  A. 

FINITE  AMPLITUDE  ACOUSTICS  AND  THE  1812  OVERTURE 

There  are  few  outdoor  musical  events  which  compare  to  a  live  performance  of 
Peter  Ilyich  Tchaikovsky’s  Overture  1812  with  cannons.  This  Opus  No.  49,  universal 
in  its  appeal,  has  packed  audiences  for  orchestras  J1  over  the  world  since  its  first 
performance  in  1882.  The  score  calls  for  not  only  cannons,  but  also  the  pealing  of 
church  bells  and  a  brass  band  ad  libitum. 

This  appendix  will  explore  some  of  the  aspects  of  1812  performances  with  can¬ 
nons  as  they  relate  to  the  subject  matter  of  this  dissertation,  the  propagation  of 
finite  amplitude  waves  outdoors.  For  example,  are  the  cannon  sounds  loud  enough 
to  exhibit  finite  amplitude  behavior?  Could  the  amplitudes  produced  be  high 
enough  to  change  the  propagation  time  of  the  cannon  blast  sounds,  interfering  with 
a  clean  performance  of  the  music?  Or  could  the  peak  pressures  be  so  loud  as  to 
cause  pain  for  certain  individuals? 

This  work  is  420  measures  long  and  runs  approximately  15  minutes,  if  no  cuts 
are  made  to  the  score.  Sixteen  cannon  blasts  were  specified  by  the  composer  in  the 
piece,  5  shots  in  measures  328  to  332,  and  11  from  measures  388  to  404.  Conduc¬ 
tors  often  take  substantial  liberties  with  the  composer’s  score,  however.  The  author 
has  heard  some  performances  where  a  single  shot  is  put  in  on  the  first  beat  of  meas¬ 
ure  36.  Other  performances  omit  Tchaikovsky’s  first  group  of  5  shots  entirely  due 
to  their  difficult  syncopated  rhythm,  and/or  continue  to  blast  away  between  meas¬ 
ures  405  and  the  end.131 
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Getting,  the  blasts  right  on  the  beat  in  time  with  the  orchestra  is  a  major 
difficulty,  as  one  might  expect.  A  chain  of  command  between  a  conductor  and 
some  national  guardsmen  firing  blanks  from  howitzers  probably  will  not  get  the  job 
done.  The  composer  originally  hoped  that  electrical  ignition  could  be  used,  a  con¬ 
ductor  having  16  switches  to  flip  at  the  appropriate  moments.132  A  technique  similar 
to  this  has  been  used  for  several  years  now  by  J.  Paul  Barnett  of  South  Bend  Repli¬ 
cas,  Inc.  Mr  Barnett,  a  builder  of  Arable  antique  cannons  for  museums  and  historic 
landmarks,  has  developed  an  electronic  system  where  he  will  fire  16  Lyle  guns  for 
an  orchestra  safely  and  on  the  beat  These  guns  originally  were  used  to  shoot  life 
saving  lines  from  shore  to  ships  in  distress. 

The  author  was  allowei  to  observe  and  record  one  of  Mr.  Barnett's  12  or  so 
concerts  per  summer,  a  performance  with  the  Chicago  Symphony  Orchestra  at  the 
Ravinia  Summer  Music  Festival  in  Highland  Park,  Illinois  on  August  13,  1989.  A 
0.0064  m  Larson  Davis  microphone  was  mounted  0.91  m  abr  the  ground,  12.2  m 
from  the  closest  gun  and  16.2  m  from  the  furthest  gun.  The  guns  were  in  two 
rows,  all  facing  away  from  the  recording  apparatus.  The  microphone  signal  was 
preamplified  and  fed  into  an  Nagra  DJ  reel-to-reel  recorder  at  0.38  m/s  speed.  The 
recording  apparatus  was  calibrated  with  a  Bruel  and  Kjaer  Pistonphone  with  a  refer¬ 
ence  level  of  123.9  dB  sound  pressure  level.  Mr.  Barnett  used  0.092  kg  black 
powder  charges,  wrapped  in  aluminum  foil  and  plastic  wrap.  There  was  a  thunder¬ 
storm  just  prior  the  the  concert,  and  the  grassy  ground  between  the  recording  dev¬ 
ices  and  the  guns  was  water  soaked  during  the  performance. 
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In  the  laboratory  the  blasts  were  analyzed  using  a  Hewlett  Packard  1727A 
Storage  Oscilloscope.  Based  on  the  voltage  waveforms  from  the  blasts,  referenced 
to  the  pistonphone  source,  the  following  peak  sound  pressure  levels  in  dB  were 
found  for  the  16  shots  in  the  order  they  were  fired:  141.8,  138.5,  136.5,  137.9, 
139.6,  135.0,  138.5,  139.0,  136.9,  140.5,  139.0,  139.8,  138.5,  138.5,  141.0,  and 
141.0.  In  most  of  the  blasts  two  peaks  were  observed,  an  initial  large  direct  wave 
peak  (SPL  just  given)  and  a  smaller  peak  occurring  later  due  to  some  reflection. 
The  blast  wave  forms  including  reflection  were  24  ms  to  32  ms  in  length. 

With  these  data  in  hand,  it  is  possible  to  speculate  on  whether  finite  amplitude 
propagation  effects  have  a  significant  impact  on  these  1812  pulses.  Considering  that 
the  sound  was  roughly  falling  off  as  1/r,  it  is  clear  that  finite  amplitude  nonlinear 
effects  should  be  important  only  in  the  near  field  of  the  guns.  If  the  level  has  al¬ 
ready  fallen  to  a  141  dB  maximum  at  a  minimum  distance  of  12.192  m,  no 
significant  change  in  the  speed  of  propagation  would  be  noticed.  Adding  a  few  dB 
to  account  for  the  receiving  equipment  being  behind  the  gun  would  make  little 
difference.  One  can  predict  that  using  slightly  more  black  powder  for  larger  charges 
would  not  change  this  hypothesis. 

On  the  other  hand,  since  many  individuals’  threshold  of  pain  is  below  140  dB, 
the  peak  pressures  from  Barnett’s  Lyle  guns  could  be  uncomfortably  loud  even 
several  hundred  feet  from  them.  It  is  important  to  remember  that  in  a  large  crowd 
of  people  with  no  ear  protection,  some  individuals  will  be  more  sensitive  to  the  blast 
noise  than  others,  and  that  the  audience  and  musicians  should  be  kept  a 
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respectable  distance  from  these  guns.  At  Ravinia,  some  audience  members  were 
approximately  30.5  m  from  the  closest  gun,  which  is  probably  too  close. 
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APPENDIX  B. 

COMPUTER  PROGRAM  DESCRIPTION 

This  appendix  describes  the  computer  program,  mcnalp,  employed  in  this 
thesis.  The  acronym  mcnalp  stands  for  MacCormack  Nonlinear  Air-Linear  Porous 
and  uses  the  second-order  in  time  and  fourth-order  in  space  MacCormack  finite 
difference  schemes  for  the  main  part  of  the  nonlinear  air  equations  in  two- 
dimensions  and  the  second-order  in  time  and  space  MacCormack  method  for  the 
linear  pore  equations  in  one-dimension.  The  code  is  over  5000  lines  of  FOR- 
TRAN77  statements  and  was  optimized  for  the  Cray  type  supercomputers.  The 
language  FORTRAN  was  used  since  that  was  the  only  scientific  language  available 
on  the  Cray  other  than  assembler  in  the  beginning  of  this  research. 

The  program  may  be  run  in  its  current  form  to  model  blast  waves  propagating 
outdoors  either  in  the  free  field  or  over  a  hard  or  a  porous  model  of  the  ground. 
One  can  also  make  minor  modifications  to  the  program  to  model  the  propagation  of 
electric  spark  pulses  in  air  or  weak  explosions  under  water. 

As  it  is  described  here,  the  program  is  ready  for  running  in  a  batch  mode  Cray 
environment  All  of  the  subroutines  are  in  one  file  except  for  one  routine 
‘'d8pimg”  which  physically  writes  out  a  National  Center  for  Supercomputing  Appli¬ 
cations  (NCSA)  Hierarchical  Data  Format  (HDF)  raster  image  file.  The  HDF  rou¬ 
tines  are  public  domain  software  routines  available  from  NCSA..  To  change  various 
parameters  which  control  the  operation  of  mcnalp,  the  user  must  modify  the  code 
itself  in  the  main  program.  At  a  later  date,  this  drawback  will  be  eliminated  by 


148 


providing  an  interactive,  user  friendly  front  end.  For  this  description,  however,  the 
current  batch  mode  implementation  will  be  assumed. 

The  program  mcnalp  provides  several  types  of  output,  some  of  which  are  more 
useful  than  others  for  particular  applications.  The  program  produces  image  files  of 
raster  images  in  HDF  form,  numerical  receiver  files  of  pressure  versus  time  at 
specified  locations,  and  a  file  of  maximum  pressures  encountered  along  the  ground 
surface  called  rpmx.out  These  are  the  three  most  often  used  types  of  output 
mcnalp  provides.  The  program  can  also  produce  ASCII  output  files  of  pressures  if 
the  raster  data  output  is  not  desired  and  can  print  almost  all  of  the  numbers  com¬ 
puted  in  the  program  using  a  debugging  option.  This  last  type  of  output  was  used 
only  in  the  early  stages  of  this  research. 

Nine  logical  variables,  i.e.,  true  or  false  variables,  control  the  major  ways  the 
program  can  operate.  When  “debug”  is  true,  debugging  output  is  printed  to  the 
screen.  If  “raster”  is  true,  the  HDF  image  files  are  produced.  Otherwise,  ASCII 
output  is  placed  in  the  file  ascii.out  If  “sscale”  is  true  the  HDF  files  are  printed 
with  the  variables  scaled  by  the  factor  1/r  which  keeps  the  images  from  fading  away 
for  long  distance  propagation.  The  flag  “allrasf  ’  specifies  whether  all  the  variables 
will  be  printed  in  raster  HDF  files,  or  just  the  acoustic  pressures.  The  variable 
“rcvrs”  indicates  whether  the  numerical  receiver  data  will  be  produced.  The  loca¬ 
tion  of  the  receivers  must  be  specified  by  the  user  in  moialp^s  main  subprogram. 
The  flag  “rowpmax”  tells  mcnalp  whether  or  not  to  save  the  values  of  maximum 
pressure  recorded  along  a  line  in  the  d  direction,  usually  the  ground  surface.  If  this 
flag  is  turned  on,  the  user  must  specify  the  distance  between  the  source  and  this  line 
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of  maximum  pressure  data  in  the  main  subprogram.  The  two  variables  “gnd”  and 
“hard”  together  tell  whether  a  ground  surface  is  used  and  if  so  what  kind  is  used. 
If  “gnd”  is  false,  a  free  field  calculation  is  made.  If  “gnd”  is  true,  “hard”  being 
true  specifies  that  a  hard  ground  is  in  use,  and  “hard”  being  false  gives  the  porous 
ground  case.  The  last  variable  “doupdat”  indicates  whether  updating  in  the  d  direc¬ 
tion  is  used  or  not.  If  “doupdat”  is  true,  an  updating  routine  is  called  every  few 
time  steps  to  conserve  computer  memory  by  overwriting  a  part  of  the  memory  that 
the  wave  leaves  by  more  free  field  that  the  wave  can  propagate  into. 

The  coordinate  system  used  in  this  program  is  cylindrical.  The  program  em¬ 
ploys  some  different  names  for  the  coordinate  axes  compared  to  the  rest  of  this 
thesis,  however.  While  in  the  main  text  of  this  thesis  the  coordinate  system  used 
was  {d,  z,  0)  with  0  variations  eliminated,  mcnalp  currently  uses  the  two  notations 
(r,  z ,  0)  and  (y  1,  y 2,  0).  Therefore,  whenever  the  reader  sees  ay  1  or  a  cylindri¬ 
cal  r  in  the  program,  it  always  means  d.  Similarly  y  2  always  means  z . 

The  main  grid  array  “grid”  in  the  air  is  rectangular  in  d  and  z.  The  coordinate 
d  goes  from  y  1  =  1  to  y  1  =  numyls,  and  z  goes  from  y2  =  1  to  y2  =  numy2s. 
The  array  “grid”  contains  the  variables  p\  u,  u,  and  Sfr  at  the  present  time  n  and 
the  temporary  time  h  used  in  the  MacCormack  finite  difference  methods.  The  vari¬ 
ables  p  and  T'  are  kept  in  arrays  of  the  same  names.  The  porous  medium  is  kept  in 
an  array  “pore”  which  contains  both  p'  and  v  for  the  pores  at  the  times  n  and  h. 
All  of  these  variables  are  kept  in  a  Fortran  common  block  /gridvars/.  Because  the 
fourth-order  in  space  MacCormack  schemes  of  the  air  cannot  be  used  near  the 
boundaries,  /gridvars/  also  contains  the  arrays  “specr”  and  “specz”  to  help  with 
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data  storage  in  implementing  the  second-order  in  space  MacCormack  schemes  near 
the  boundaries. 

The  other  common  blocks  in  this  program  hold  global  variables  grouped  togeth¬ 
er  for  convenience.  The  block  /  times/  contains  information  needed  for  printing  the 
raster  images  at  specific  times,  the  values  of  and  the  time  at  which  the  program 
will  stop.  The  block  /cons/  holds  the  constants  needed  for  air  propagation  such  as 
the  sound  speed  and  the  other  thermodynamic  variables.  On  the  other  hand, 
/peons/  contains  the  constants  needed  for  the  porous  medium.  The  massive  block 
/bndvars/  groups  together  the  variables  needed  for  the  boundary  calculations.  Since 
“grid”  is  rectangular,  four  boundaries  are  needed,  the  boundary  close  to  the  source, 
the  outer  boundary  opposite  the  source  boundary,  the  ground  boundary,  and  the 
upper  boundary. 

Some  smaller  common  blocks  are  much  more  limited  in  their  scope.  The  block 
/flags/  holds  the  logical  variables  described  above.  The  block  /diffcon/  contains  the 
constants  used  in  the  main  MacCormack  routines  which  are  precomputed  to  speed 
program  execution.  The  common  block  /pulsvars/  holds  the  variables  to  describe 
the  modified  Jack  Reed  pulse  for  the  initial  conditions  of  the  simulation.  The  infor¬ 
mation  to  control  the  updating  in  the  d  direction  routines  is  contained  in  the  block 
/upblk/.  The  block  /rcvr/  holds  numerical  receiver  height  and  distance  informa¬ 
tion.  Finally,  /rpmax/  includes  the  variables  to  save  the  maximum  pressures  along 
a  line  in  the  d  direction,  usually  along  the  ground  surface. 

The  subroutines  in  the  program  are  for  the  most  part  self-documenting  with 
ample  comment  lines.  The  largest  routines  “diffarp,”  ‘‘diffarc,”  “diffazp,” 
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“diffazc,”  “diffbrp,”  “difibrc,”  “diffbzp,”  and  “diffbzc”  were  the  most  time  con¬ 
suming  to  write  and  perform  the  MacCormack  time  difference  methods  in  the  air. 
The  diffa  routines  use  the  “a”  version  of  the  MacCormack  methods,  while  the  diffb 
routines  use  the  “b”  version.  The  next  to  last  of  the  letters  in  the  names  of  these 
eight  routines  specifies  first  what  direction  the  routine  is  for,  either  r,  i.e.,  d,  or  z, 
and  the  last  letter  tells  whether  this  routine  gives  the  predictor  step  or  the  corrector 
step.  For  example,  “diffbzp”  is  the  predictor  step  of  the  b  MacCormack  method  in 


the  z  direction. 


The  boundary  conditions  are  implemented  in  routines  “upbndl,”  “gndbndl,” 
“outbndl,”  “srcbndl,”  “prebnd,”  and  “bnd2.”  The  first  four  routines  compute 
and  store  the  information  needed  for  the  boundaries  under  free  field  conditions, 
and  “bnd2”  applies  either  these  data  and/or  the  appropriate  ground  condition. 
Routine  “prebnd”  precalculates  values  needed  for  the  numerical  boundary  condi¬ 
tions.  Subroutine  “porous”  computes  all  of  the  variable  values  in  the  porous  medi¬ 
um  for  use  as  a  porous  boundary  condition  in  “bnd2.”  The  equations  of  state  in  the 
air  at  the  times  n  and  n  are  used  in  the  routines  “pandt”  and  “pandt2,”  respective¬ 
ly,  to  give  updated  values  of  p  and  T'.  Subroutines  “enta”  and  “entb”  are  used  in 
addition  to  the  “diffXXX”  modules  to  keep  the  values  of  Sfr  correctly  updated  near 
the  boundaries  in  the  air. 


Two  of  the  unique  features  of  mcnalp  include  saving  memory  by  updating  or 
windowing  in  the  d  direction  and  recording  the  maximum  pressure  values  along  the 
ground.  The  routine  “update”  performs  the  first  function,  and  “getrpmax”  the 
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second.  The  subroutine  “  wrpmax”  writes  out  the  file  rpmx  out  if  the  maximum 
pressure  information  is  desired. 

The  other  routines  which  give  results  for  later  analysis  are  “pmtifok,” 
“pmtgrid,”  and  “prraster.”  The  program  calls  ‘‘pmtifok”  every  time  step  to  see  if 
any  images  should  be  saved  and  to  save  the  numerical  receiver  waveforms,  should 
there  be  any  receivers  specified.  If  the  raster  data  output  is  needed,  “pmtifok”  calls 
“prraster,”  or  if  ASCII  images  are  needed,  it  calls  “pmtgrid.” 

The  first  of  the  rest  of  the  Fortran  modules  called  by  the  main  program  are 
“predifF’  which  precomputes  the  constants  needed  for  the  MacCormack  routines 
“diffXXX.”  The  routines  “getw,”  “getw2,”  “getwt,”  and  “getwt2”  compute  one- 
directional  derivatives  so  the  “diffXXX"  routines  can  calculate  the  mixed  deriva¬ 
tives  which  occur  in  the  dissipative  terms  of  the  air  equations.  This  is  one  place 
where  one  might  be  able  to  improve  mcnalp,  by  incorporating  these  routines  into 
the  “diffXXX”  routines.  The  program  also  calls  the  real  functions  “rm2”  and 
“rpot2”  which  are  the  modified  Reed  blast  pulse  pressure  and  potential  waveforms, 
respectively. 

The  main  routine  in  mcnalp  can  be  thought  of  in  four  sections  after  all  the 
variable  declarations.  In  the  first  section  the  user  specifies  the  problem  to  be  solved 
by  setting  the  flags,  defining  the  constants  to  be  used,  placing  the  boundaries,  speci¬ 
fying  the  times  to  print  images,  and  placing  the  numerical  receivers.  The  second 
section  of  mcnalp  generates  the  discretized  grid  on  which  the  problem  will  be 
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solved,  and  the  third  section  initializes  the  variables  to  their  initial  condition  for  the 
beginning  of  the  calculation. 

The  last  section  of  mcnalp  is  the  main  loop  of  the  program  in  which  the  nu¬ 
merical  solution  is  advanced  in  time.  For  every  iteration  of  this  loop  the  variables 
are  all  stepped  in  time  one  Jtf.  Each  time  the  loop  is  executed,  two  paths  are  alter¬ 
nately  taken.  This  alternation  is  necessary  since  the  McCormack  “a’;  and  “b” 
schemes  must  be  alternated,  as  well  as  the  r,  i.e.,  d ,  and  z  directions. 

The  first  path  through  the  loop  is  now  described.  The  second  path  is  similar. 
The  program  mcnalp  first  computes  the  radiation  boundary  conditions  and  saves 
these  numbers  using  “upbndl,”  “gndbndl,”  “outbndl,”  and  “srcbndl.”  The  r, 
i.e.,  d ,  direction  routines  are  then  called.  The  routine  “getwt”  precomputes  infor¬ 
mation  for  “diffarp”  which  comes  next  Subroutine  “enta”  then  updates  the  entro¬ 
py  variable  Sfr  near  the  boundaries.  Routines  “pandt2”  and  “getwt2”  then  update 
p  and  T'  at  the  temporary  time  and  precompute  information  at  the  temporary  time 
for  “diffarc”  which  follows.  An  invocation  of  to  obtain  the  pressure  and 

temperature  deviation  completes  the  r  direction  par:  .J  the  main  loop’s  first  path. 

The  z  direction  part  of  the  main  loop’s  first  path  begins  by  “getw”  precomput¬ 
ing  information  for  “diffazp”  which  follows  it  The  routine  “pandt2”  then  updates 
p  and  T'  *  the  temporary  time.  Subroutine  “getw2”  then  precomputes  values  for 
“diffazc”  which  comes  next  The  “bnd2”  routine  is  then  called  to  update  the 
boundary  values  for  the  air  grid.  If  the  porous  ground  is  employed  in  the  calcula¬ 
tion,  “bnd2”  calls  the  routine  “porous”  at  this  time.  After  “bnd2”  finishes,  the 
routine  “pandt”  updates  p  and  T'  at  this  time.  The  program  has  now  advanced  the 
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numerical  solution  one  time  step,  Jtf.  This  pass  of  the  main  loop  ends  by  calling 
“getrpmax,”  “pmtifok,”  and  “update.” 

In  the  next  iteration  of  the  main  loop  the  “b”  routines  are  called,  first  in  the  z 
direction  then  in  the  r,  i.e.,  d,  direction.  It  then  also  calls  “getrpmax,”  “pmtifok,” 
and  “update.”  Next  the  “a”  path  of  the  main  loop  is  executed  again.  This  process 
continues,  the  program  alternately  passing  through  the  two  paths  of  the  main  loop, 
until  mcnalp  terminates  when  the  time  reached  exceeds  the  maximum  time  the 
user  has  specified. 

The  last  mqjor  component  of  mcnalp  which  requires  detailed  description  is  the 
updating  procedure.  If  the  flag  “doupdat”  is  false,  no  updating  is  performed,  and 
the  grid  variables  span  the  d  direction  from  “minyl”  to  “miny2”  using  the  discreti¬ 
zation  i  going  from  1  to  “numyls.”  The  problem  here  is  that  a  wave,  initially 
between  “minyl”  and  “maxyl,”  moving  in  the  d  direction  will  eventually  pass  the 
“maxyl”  boundary  and  propagate  out  of  the  numerical  domain. 

If  the  flag  “doupdat”  is  true,  however,  a  wave  moving  in  the  d  direction  will 
never  move  out  of  the  numerical  domain  because  every  few  time  steps  the  domain 
is  extended  in  the  +d  direction  while  simultaneously  deleted  in  the  —d  direction. 
Performing  these  extensions  beyond  the  outer  boundary  and  deletions  close  to  the 
source  boundary  windows  a  pulse  moving  in  the  Hi  direction. 

This  action  saves  much  computer  memory  for  a  long  distance  propagation  run. 
Instead  of  having  all  of  the  distance  between  the  initial  pulse  position  and  its  final 
position  in  memory  at  once,  only  the  window  that  the  pulse  is  in  at  any  one  time  is 
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kept  in  memory. 

The  windowing,  or  updating,  procedure  is  performed  with  the  array  “ii.”  In 
mcnalp,  whenever  a  variable  is  given  d  or  z  positions  such  as  p(i,j)  for  the  pressure 
at  d  =  “yl(i)”  and  2  =  “y2(j),”  the  notation  p(ii(i),j)  is  used.  If  “doupdat”  is 
false,  no  updating  is  performed  and  the  “ii”  array  takes  the  values  “ii(l)”=  1, 
“ii(2)”  =  2,  .  .  .,  “ii(n)”  =  n,  etc.  Thus  in  the  “doupdat”  false  case  the  array  “ii” 
acts  as  if  it  is  not  there  at  all. 


If  “doupdat”  is  true,  however,  after  the  wave  has  propagated  a  certain  dis¬ 
tance,  an  update  is  performed.  Suppose  there  were  ten  points  in  the  d  direction, 
“numyls”=  10,  and  “ii”  and  “yl”  had  the  values  “ii(l)”=  1,  “ii(2)”=  2,  .  .  ., 
“ii(  10)”=  10  and  “yl(ii(l))”=  1.0  m,  “yl(ii(2))”=  1.1  m,  “yl(ii(3))”=  1.2  m,  .  . 
.,  “yl(ii(10))”=  1.9  m,  respectively,  initially.  If  the  update  shifted  things  by  2 
points,  one  would  be  left  with  “ii(l)”=  3,  “ii(2)”=  4,  .  .  .,  “ii(8)”=  10, 


“ii(9)”=  1,  “ii(10)”=  2  and  “yl(ii(l))”=  1.2  m,  “yl(ii(2))”=  1.3  m,  .  .  ., 


“yl(ii(8))”=  2.0  m,  “yl(u(9))”=  2.1  m,  “yl(ii(  10))”=  2.2  m.  Here  the  d  dis¬ 
tances  1.0  m  and  1.1  m  are  no  longer  represented,  having  been  replaced  by  dis¬ 


tances  2.0  and  2.1  m.  By  repeating  this  method  over  and  over  for  updates  as  the 
wave  moves,  the  pulse  is  windowed  for  propagation  in  the  d  direction. 


To  implement  this  updating  procedure,  the  user  needs  to  specify  a  variable  in 
the  main  program,  “ng.”  The  variable  “ng”  is  the  number  of  grid  increments,  i.e., 
the  number  of  Ax’s,  the  wave  should  move  before  an  update  is  performed.  Anoth¬ 
er  variable  the  program  automatically  initializes  for  the  updating  procedure  is  “ut.” 
The  variable  “ut”  is  the  update  time  and  gives  the  amount  of  time  elapsed  since  the 
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last  update.  The  test  “Is  ut  c0  >  ng  Ar  ?”  is  made  at  the  bottom  of  the  main  loop 
of  the  program,  and  routine  “update’’  is  called  if  this  test  is  true.  The  routine  “up¬ 
date”  zeros  the  variable  “ut?’  every  time  the  procedure  is  invoked  to  prepare  for 
the  next  update.  The  subroutine  “ut”  is  incremented  by  on  every  iteration  of 
the  main  loop  to  count  up  until  the  update  test  is  true. 

One  final  note  about  the  present  version  of  mcnalp.  Statements  reading 
CDIR$  IVDEP  were  placed  before  some  of  the  do  loops  in  the  subroutines.  Under 
Cray  FORTRAN  the  statement  is  a  compiler  directive  to  Ignore  Vector  DEPenden- 
cies.  These  compiler  directives  were  inserted  to  tell  the  Cray  FORTRAN  compiler 
that  even  though  the  do  loop  might  look  as  though  it  could  contain  a  dependency,  it 
does  not  and  should  be  vectorized.  Dependencies  are  sequences  of  memory  refer¬ 
ence  operations  that  keep  a  loop  from  vectorizing,  and  eliminating  dependencies 
gives  a  marked  increase  in  performance.  Under  most  FORTRAN  compilers  these 
CDIR$  IVDEP  statements  are  interpreted  as  comment  statements. 

To  give  some  idea  of  how  well  mcnalp  vectorizes  and  performs  on  the  Cray 
series,  a  short  run  was  made  on  both  the  Cray-2  at  NCSA  and  the  CERL  Pyramid 
Pyr90x  mainframe  computer.  Global  optimization  was  used  with  the  Fortran  com¬ 
pilers  on  both  systems.  On  the  CERL  Pyramid  with  no  one  else  logged  in,  the  pro¬ 
gram  took  approximately  90  min  of  CPU  time.  The  Cray-2  ran  the  same  program 
in  5.3  s,  approximately  a  speed-up  of  1000  times. 
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